$IBJOB         MAP
$IBFTC ZERO                                                                  000
      BLOCK DATA                                                             001
      COMMON/ERR/I1/NAMEA/I2(10)/NAME1/I3(5)/NAMEW/I4(3)/NAME2/I5(293)       002
     X/NAME3/I6(99)/NAME4/I7(50)                                             003
      DATA I1/0/,(I2(I),I=1,10)/0,0,0,0,0, 0,0,0,0,0/,                       004
     X(I3(I),I=1,5)/0,0,0,0,0/,(I4(I),I=1,3)/0,0,0/,(I5(I),I=1,293)/         005
     X0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0,      006
     X0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0,      007
     X0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0,      008
     X0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0,      009
     X0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0,      010
     X0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0,      011
     X0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0,      012
     X0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0,      013
     X0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0,      014
     X0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0/                     015
      DATA (I6(I),I=1,99)/                                                   016
     X0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0,      017
     X0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0,      018
     X0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0,      019
     X0,0,0,0,0, 0,0,0,0/                                                    020
      DATA (I7(I),I=1,50)/                                                   021
     X0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0,      022
     X0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0/                            023
      END                                                                    024
$IBFTC N1MAIN                                                           N1MA0000
      DIMENSION   A(10)    ,BS(10)   ,WPWG(10)           ,PP(50,100),   N1MA0001
     X PS(50,50) ,PTB(50,8),PW(150)  ,R(10)    ,TRL1(5)  ,TRL2(5)   ,   N1MA0002
     X TRL3(5)   ,TRL4(5)  ,TRM1(5)  ,TRM2(5)  ,TRM3(5)  ,TRM4(5)   ,   N1MA0003
     X FA(100)   ,WDP(10)  ,GF(10)                                      N1MA0004
     X,IT(20)    ,KT(20)                                                N1MA0005
     X,CK(10)    ,CL(10)   ,CKP(10)  ,WPWT(10) ,R234(10) ,XXX(31)       N1MA0006
     X,SAUR(6)                                                          N1MA0007
      COMMON CAPN  ,CPL   ,CPS   ,CTOP  ,DL    ,DR    ,INSERT,DTWI  ,   N1MA0008
     XEW    ,WPWG  ,GAMMA ,GMG0  ,HPL   ,HPS   ,IMAX  ,N1    ,N2    ,   N1MA0009
     XPS    ,PG0   ,PW    ,R     ,RCAP  ,RHOP0 ,SMP   ,TG0   ,TPM   ,   N1MA0010
     XA     ,BS    ,CC    ,CONST ,EOW1  ,EOW2  ,HIT   ,ICNT  ,IPP   ,   N1MA0011
     XJMAX  ,MW           ,PP    ,PTB   ,TRL1  ,TRL2  ,TRL3  ,TRL4  ,   N1MA0012
     XTRM1  ,TRM2  ,TRM3  ,TRM4  ,LRC   ,DELTA ,FA    ,VGR   ,CPG   ,   N1MA0013
     XNILP  ,THI   ,GF    ,WDG   ,IFG   ,RHOG0 ,THRUST,TI    ,WDP       N1MA0014
      COMMON TMD   ,IPPX  ,BC    ,                                      N1MA0015
     XTHII  ,DZMIN ,ITOT  ,ZAX   ,RRT   ,RT    ,VK    ,T     ,CG    ,   N1MA0016
     XA0C   ,DZI   ,ZI    ,ZJ    ,NZONES,CK    ,CL    ,CKP   ,R1    ,   N1MA0017
     XZ1    ,UG    ,VG    ,TT0   ,DG    ,UGR   ,UGZ   ,VRG   ,VGZ   ,   N1MA0018
     XZW    ,R234  ,WPWGT ,WPWT  ,XXX   ,THJ   ,SAUR  ,EPS   ,RWMAX ,   N1MA0019
     XZWMAX ,THIW  ,THJW  ,XXXXX ,THFD  ,THJD  ,IP    ,IT    ,KT    ,   N1MA0020
     XLAST  ,ZMAX                                                       N1MA0021
C                                                                       N1MA0022
        COMMON/ERR/IERR                                                 N1MA0023
      COMMON/NAMEA/VAR(6),IORDER,SOWCON,ISW1,ISW2                       N1MA0024
      COMMON/NAME1/THF,MARK,UUS,UDUDS,ZSINK                             N1MA0025
      COMMON/NAMEW/IWALL,THIR,THJR                                      N1MA0026
C                                                                       N1MA0027
      DIMENSION FJ(59)                                                  N1MA0028
     *,IA(8460)
      EQUIVALENCE (IA(1),CAPN)
      DATA(FJ(I),I=1,59)/0.,0.,1.,1.25,1.,1.255,1.,1.26,1.001,1.265     N1MA0029
     X ,1.002,1.582                                                     N1MA0030
     1 ,1.063,1.995,1.141,2.51,1.224,3.16,1.315,3.98,1.412,5.01,1.517   N1MA0031
     2 ,6.31,1.625,7.95,1.745,10.,1.874,12.6,2.026,15.82,2.186,19.95    N1MA0032
     3 ,2.364,25.1,2.555,31.6,2.76,39.8,3.,50.1,3.252,63.1,3.534        N1MA0033
     4 ,79.5,3.825,100.,4.155,316.,7.9,1000.,20.,1001.,20.02,1.E5       N1MA0034
     5 ,2000.,0.,0./                                                    N1MA0035
C                     1   2   3    4   5  6   7    8   9                N1MA0036
      NAMELIST/DATA/CAPN,CPL,CPS,PR  ,DL,DR      ,DTWI,EW,              N1MA0037
C     10   11   12   13  14  15  16 17 18  19 20                        N1MA0038
     XFA,GAMMA,GMG0,HPL,HPS,IMAX,N1,N2,PS,PC ,PW,                       N1MA0039
C     21  22   23   24  25  26  27  28   29   30  31                    N1MA0040
     XR,RCAP,RHOP0,SMP,TG0,TPM,CPG,NILP,THID,WPWT,RT,                   N1MA0041
C      32   33    34   35  36  37 38 39  40                             N1MA0042
     XRRT,WPWGT,DZMIN,IFG,TMD,ZAX,ZI,ZJ,DZI,                            N1MA0043
C       41     42    43   44   45   46  47   48   49                    N1MA0044
     XIORDER,SOWCON,ISW1,ISW2,ISW3,VAR,FLAG,THJD,SAUR,                  N1MA0045
C      50   51   52    53    54   55   56                               N1MA0046
     XTHFD,EPS,RWMAX,ZWMAX,IWALL,THIW,THJW                              N1MA0047
C                                                                       N1MA0048
      DO 20 IB=1,8460
      IA(IB)=0
 20   CONTINUE
      DO 10 I=1,59                                                      N1MA0049
 10   FA(I)=FJ(I)                                                       N1MA0050
      DL=.2                                                             N1MA0051
      DR=.2                                                             N1MA0052
      DTWI=3.                                                           N1MA0053
      DZI=.002                                                          N1MA0054
      DZMIN=.002                                                        N1MA0055
      EW=.001                                                           N1MA0056
      IMAX=5                                                            N1MA0057
      IORDER=3                                                          N1MA0058
      ISW1=0                                                            N1MA0059
      ISW2=0                                                            N1MA0060
      NILP=15                                                           N1MA0061
      N1=1000                                                           N1MA0062
      N2=1                                                              N1MA0063
      SAUR(1)=-.15                                                      N1MA0064
      SAUR(2)=1.                                                        N1MA0065
      SAUR(3)=.5                                                        N1MA0066
      SAUR(4)=.3                                                        N1MA0067
      SAUR(5)=-1.                                                       N1MA0068
      THJD=9.                                                           N1MA0069
      THFD=5.                                                           N1MA0070
      THIW=12.                                                          N1MA0071
      TMD=1000.                                                         N1MA0072
      VAR(1)=.3                                                         N1MA0073
      VAR(2)=0.                                                         N1MA0074
      VAR(3)=0.                                                         N1MA0075
      VAR(4)=.1                                                         N1MA0076
      VAR(5)=.1                                                         N1MA0077
      ZI=3.                                                             N1MA0078
      ZJ=2.                                                             N1MA0079
C                                                                       N1MA0080
 1    CONTINUE                                                          N1MA0081
      IERR=0                                                            N1MA0082
      DO 11 I=1,50                                                      N1MA0083
      DO 11 J=1,50                                                      N1MA0084
 11   PS(I,J)=0.                                                        N1MA0085
C                                                                       N1MA0086
      READ(5,DATA)                                                      N1MA0087
C                         4633.056=144.*32.174  (CONVERT FROM PSI)      N1MA0088
      PG0=PC*4633.056                                                   N1MA0089
      CTOP=CPG/PR                                                       N1MA0090
 12   DO 2 I=1,11                                                       N1MA0091
      ITOT=I-1                                                          N1MA0092
      IF(R(I)) 2,3,2                                                    N1MA0093
 2    CONTINUE                                                          N1MA0094
 3    RHOG0=PG0/RCAP/TG0                                                N1MA0095
      RHOP0=RHOG0*WPWGT                                                 N1MA0096
      DO 4 I=1,ITOT                                                     N1MA0097
 4    WPWG(I)=WPWGT*WPWT(I)                                             N1MA0098
      VK=GAMMA*(1.+WPWGT*CPL/CPG)/(1.+WPWGT*CPL/CPG*GAMMA)              N1MA0099
      BC=(1.+WPWGT)/(1.+WPWGT*CPL/CPG)                                  N1MA0100
      WDG=RHOG0*(2./(VK+1.))**(1./(VK-1.))*SQRT (2.*CPG*TG0             N1MA0101
     X    /BC*(VK-1.)/(VK+1.))*3.1415927*RT*RT                          N1MA0102
      T=1./RRT                                                          N1MA0103
      THI=.017453293*THID                                               N1MA0104
      THJ=.017453293*THJD                                               N1MA0105
      THF=.017453293*THFD                                               N1MA0106
      THIR=.017453293*THIW                                              N1MA0107
      A0C=4.5/SMP*GMG0*RT                                               N1MA0108
      CG=3.42/SQRT (GAMMA)*CTOP/CPG                                     N1MA0109
      WRITE (6,90)                                                      N1MA0110
      PW(1)=1.+RRT*(1.-COS (THIR))                                      N1MA0111
      PW(2)=RRT*SIN (THIR)                                              N1MA0112
      CALL PARTIL                                                       N1MA0113
      IF(IERR .NE. 0) GO TO 1                                           N1MA0114
C                                                                       N1MA0115
      CC=GMG0/(TG0*RCAP)**CAPN                                          N1MA0116
      CONST=4.5/SMP*CC*RT                                               N1MA0117
      CG=3.42/SQRT (GAMMA)*CTOP/CPG                                     N1MA0118
      BSC=3.*CTOP/SMP*CC*RT                                             N1MA0119
      DO 5 I=1,ITOT                                                     N1MA0120
      BS(I)=BSC/R(I)**2                                                 N1MA0121
 5    WDP(I)=WPWG(I)*WDG                                                N1MA0122
      THIR=.017453293*THIW                                              N1MA0123
      THJR=.017453293*THJW                                              N1MA0124
      DELTA=.005                                                        N1MA0125
      JMAX=0                                                            N1MA0126
      CALL N2MAIN                                                       N1MA0127
      WRITE (6,91)IP                                                    N1MA0128
      GO TO 1                                                           N1MA0129
C                                                                       N1MA0130
 90   FORMAT(1H1,46X,17HGAS-PARTICLE FLOW///)                           N1MA0131
 91   FORMAT(4H0IP=I3)                                                  N1MA0132
       END                                                              N1MA0133
$IBFTC FIND                                                             FIND0000
      SUBROUTINE  FIND(N,A,B,C)                                         FIND0001
      DIMENSION A(1)                                                    FIND0002
      FINDA=0.0                                                         FIND0003
      FIND1=0.0                                                         FIND0004
      FINDB=B                                                           FIND0005
      I=MAX1(A(1),2.0)                                                  FIND0006
  10  IF(FINDB.GT.A(I)) GO TO 30                                        FIND0007
      IF(I.EQ.2) GO TO 60                                               FIND0008
      I=I-2                                                             FIND0009
      GO TO 10                                                          FIND0010
  20  I=I+2                                                             FIND0011
  30  IF(A(I+2).EQ.0.0.AND.A(I+3).EQ.0.0) GO TO 60                      FIND0012
      IF(A(I+2).LE.FINDB) GO TO 20                                      FIND0013
      FIND5=A(I)                                                        FIND0014
      FIND4=A(I+2)-FIND5                                                FIND0015
      FIND1=(A(I+3)-A(I+1))/FIND4                                       FIND0016
      IF(N.EQ.0.OR.I.LT.4.OR.(A(I+4).EQ.0..AND.A(I+5).EQ.0.)) GO TO 50  FIND0017
      FIND2=FIND5-A(I-2)                                                FIND0018
      FIND3=A(I+4)-FIND5                                                FIND0019
      FINDA=(FIND1*(FIND2-FIND3)+A(I-1)-2.*A(I+1)+A(I+5))               FIND0020
     1     /(FIND4*(FIND2-FIND3)+FIND2**2+FIND3**2)                     FIND0021
  50  FIND5=FINDB-FIND5                                                 FIND0022
  60  C=(FINDA*(FINDB-A(I+2))+FIND1)*FIND5+A(I+1)                       FIND0023
      A(1)=I                                                            FIND0024
      RETURN                                                            FIND0025
      END                                                               FIND0026
$ORIGIN        LEVEL1,UT2
$IBFTC ABCALC                                                           ABCA0000
CABCALC                                                                 ABCA0001
      SUBROUTINE ABCALC                                                 ABCA0002
C                                                                       ABCA0003
      COMMON/NAMEA/VAR(6),IORDER,SOWCON,ISW1,ISW2                       ABCA0004
      COMMON/NAME1/THF,MARK,UUS,UDUDS,ZSINK                             ABCA0005
      COMMON/NAME2/VK,NFV,NVAR(6),X,U                                   ABCA0006
     1  ,ALPHA,BETA,GAMMA,DELTA,F(6),X2,X3,W(6)                         ABCA0007
     2  ,MF(6),D(6),AB(6,7),BLOCK(200),MAXIT,NUMIT                      ABCA0008
     3  ,X4,NSAV,XW,XW2,XW3,VKK,T2,XSTAR,Y                              ABCA0009
      DIMENSION A20(1),B21(1),C22(1),D23(1)                             ABCA0010
      COMMON/NAME3/                                                     ABCA0011
     1  A20,A21,A40,B21,B22,B40,B41,B42,B60,B61,B80,C22,C23,C40,C41,C42,ABCA0012
     2C43,C60,C61,C62,C63,C80,C81,C82,C100,C101,C120,D23,D24,D41,D42,D43ABCA0013
     3,D44,D60,D61,D62,D63,D64,D80,D81,D82,D83,D84,D100,D101,D102,D103  ABCA0014
     4 ,D120,D121,D122,D140,D141,D160,D40                               ABCA0015
     5 ,R,R2,R3,R4,R5,R6,R7,R8,R9,R10,R11,R12,R13,R14,R15               ABCA0016
     6 ,D1,D2,D3,D4,D5,D6,D7,D8,D9,D10,D11,D12,D13,D14,D15              ABCA0017
     7 ,UF,VF,DRDX1,DRDX2,DRDX3,DRDX4,DRDX5,DUDX1                       ABCA0018
     8 ,DUDX2,DUDX3,DUDX4,DVDX1,DVDX2,DVDX3,DVDX4                       ABCA0019
      COMMON/NAME4/                                                     ABCA0020
     *       CC0,CC1,CA201,CA102,CA120,CA220,CA320,CA121,CA321          ABCA0021
     1 ,CA122,CA140,CA340,CA141,CA160,CB202,CB103,CB220                 ABCA0022
     2 ,CB320,CB121,CB321,CB421,CB122,CB221,CB222,CB422                 ABCA0023
     3 ,CB322,CB123,CB140,CB240,CB340,CB440,CB141,CB241,CB341           ABCA0024
     4 ,CB441,CB142,CB442,CB143,CB160,CB260,CB360,CB460                 ABCA0025
     5 ,CB161,CB461,CB162,CB180,CB480,CB181,CB1100,CB120                ABCA0026
C                                                                       ABCA0027
C     FIRST ORDER EQUATION                                              ABCA0028
C                                                                       ABCA0029
      CC0=U*U-1.0                                                       ABCA0030
      CC1=1.0-VKK*U*U                                                   ABCA0031
      A20=CC0*ALPHA/(4.0*CC1)                                           ABCA0032
      A21=U*ALPHA**2/(2.0*CC1)                                          ABCA0033
      A40=U*A21*ALPHA/(8.0*CC1)                                         ABCA0034
C                                                                       ABCA0035
C     SECOND ORDER EQUATION                                             ABCA0036
C                                                                       ABCA0037
      CA201=2.0*VKK*U*ALPHA                                             ABCA0038
      B21=(4.0*CA201*A20+CC0*BETA)/(4.0*CC1)                            ABCA0039
      CA102=ALPHA**2+U*BETA                                             ABCA0040
      B22=(2.0*U*ALPHA*BETA+CA102*ALPHA+4.0*CA201*A21)/(4.0*CC1)        ABCA0041
      CA120=2.0*U*B21                                                   ABCA0042
      CA220=2.0*VKK*U*A21                                               ABCA0043
      CA320=8.0*U*A20/(VK+1.0)                                          ABCA0044
      B40=(2.0*CC0*B22+CA120*ALPHA+4.0*CA220*A20+2.0*CA320*A21)/(16.0*CCABCA0045
     11)                                                                ABCA0046
      CA121=2.0*A21*ALPHA+4.0*U*B22                                     ABCA0047
      CA321=8.0*U*A21/(VK+1.0)                                          ABCA0048
      B41=(2.0*U*A21*BETA+4.0*U*B22*ALPHA+CA121*ALPHA+4.0*CA220*A21     ABCA0049
     1      +16.0*CA201*A40+2.0*CA321*A21)/(16.0*CC1)                   ABCA0050
      CA140=A21**2+2.0*U*B41                                            ABCA0051
      CA340=16.0*U*A40/(VK+1.0)                                         ABCA0052
      B60=(4.0*U*A21*B22+CA140*ALPHA+16.0*CA220*A40+2.0*CA340*A21)      ABCA0053
     1      /(36.0*CC1)                                                 ABCA0054
      RETURN                                                            ABCA0055
      END                                                               ABCA0056
$IBFTC CCALC                                                            CCAL0000
      SUBROUTINE CCALC                                                  CCAL0001
CCCALC                                                                  CCAL0002
C                                                                       CCAL0003
      COMMON/NAMEA/VAR(6),IORDER,SOWCON,ISW1,ISW2                       CCAL0004
      COMMON/NAME1/THF,MARK,UUS,UDUDS,ZSINK                             CCAL0005
      COMMON/NAME2/VK,NFV,NVAR(6),X,U                                   CCAL0006
     1  ,ALPHA,BETA,GAMMA,DELTA,F(6),X2,X3,W(6)                         CCAL0007
     2  ,MF(6),D(6),AB(6,7),BLOCK(200),MAXIT,NUMIT                      CCAL0008
     3  ,X4,NSAV,XW,XW2,XW3,VKK,T2,XSTAR,Y                              CCAL0009
      DIMENSION A20(1),B21(1),C22(1),D23(1)                             CCAL0010
      COMMON/NAME3/                                                     CCAL0011
     1  A20,A21,A40,B21,B22,B40,B41,B42,B60,B61,B80,C22,C23,C40,C41,C42,CCAL0012
     2C43,C60,C61,C62,C63,C80,C81,C82,C100,C101,C120,D23,D24,D41,D42,D43CCAL0013
     3,D44,D60,D61,D62,D63,D64,D80,D81,D82,D83,D84,D100,D101,D102,D103  CCAL0014
     4 ,D120,D121,D122,D140,D141,D160,D40                               CCAL0015
     5 ,R,R2,R3,R4,R5,R6,R7,R8,R9,R10,R11,R12,R13,R14,R15               CCAL0016
     6 ,D1,D2,D3,D4,D5,D6,D7,D8,D9,D10,D11,D12,D13,D14,D15              CCAL0017
     7 ,UF,VF,DRDX1,DRDX2,DRDX3,DRDX4,DRDX5,DUDX1                       CCAL0018
     8 ,DUDX2,DUDX3,DUDX4,DVDX1,DVDX2,DVDX3,DVDX4                       CCAL0019
      COMMON/NAME4/                                                     CCAL0020
     *       CC0,CC1,CA201,CA102,CA120,CA220,CA320,CA121,CA321          CCAL0021
     1 ,CA122,CA140,CA340,CA141,CA160,CB202,CB103,CB220                 CCAL0022
     2 ,CB320,CB121,CB321,CB421,CB122,CB221,CB222,CB422                 CCAL0023
     3 ,CB322,CB123,CB140,CB240,CB340,CB440,CB141,CB241,CB341           CCAL0024
     4 ,CB441,CB142,CB442,CB143,CB160,CB260,CB360,CB460                 CCAL0025
     5 ,CB161,CB461,CB162,CB180,CB480,CB181,CB1100,CB120                CCAL0026
C                                                                       CCAL0027
C     THIRD ORDER EQUATION                                              CCAL0028
C                                                                       CCAL0029
      CONS1=VK+1.                                                       CCAL0030
      CONS2=VK-1.                                                       CCAL0031
      CONS3=VKK                                                         CCAL0032
C                                                                       CCAL0033
      CB202=CONS3*CA102                                                 CCAL0034
      C22=(CC0*GAMMA*.5+4.*CA201*B21+4.*CB202*A20)/(4.*CC1)             CCAL0035
      CB103=U*GAMMA/3.+ALPHA*BETA                                       CCAL0036
      C23=(4.*CA201*B22+4.*CB202*A21+CA102*BETA+CB103*ALPHA             CCAL0037
     1   +U*ALPHA*GAMMA)/(4.*CC1)                                       CCAL0038
      CB120=4.*CONS3*A20**2                                             CCAL0039
      CB220=CONS3*CA120                                                 CCAL0040
      C40=(CB120*ALPHA+4.*CB220*A20+2.*CA320*B21+2.*CC0*C22)/(16.*CC1)  CCAL0041
      CB121=4.*U*C22+2.*B21*ALPHA+8.*CONS3*A20*A21                      CCAL0042
      CB221=CONS3*CA121                                                 CCAL0043
      CB321=8.*(A20*ALPHA+U*B21)/CONS1                                  CCAL0044
      C41=(4.*CA320*B22+2.*CA321*B21+16.*CA201*B40+4.*CA220*B21         CCAL0045
     1   +CB121*ALPHA+4.*CB221*A20+4.*CB220*A21+2.*CB321*A21+CA120*BETA CCAL0046
     2   +6.*CC0*C23+4.*U*C22*ALPHA)/(16.*CC1)                          CCAL0047
      CB122=6.*U*C23+A21*BETA+4.*B22*ALPHA+4.*CONS3*A21**2              CCAL0048
      CB322=8.*(A21*ALPHA+U*B22)/CONS1                                  CCAL0049
      C42=(CB122*ALPHA+4.*CB221*A21+16.*CB202*A40+2.*CB322*A21          CCAL0050
     1   +CA121*BETA+2.*CA102*B22+16.*CA201*B41+4.*CA321*B22            CCAL0051
     2   +4.*CA220*B22                                                  CCAL0052
     3   +12.*U*C23*ALPHA+U*A21*GAMMA)/(16.*CC1)                        CCAL0053
      CB140=2.*U*C41+2.*A21*B21+16.*CONS3*A20*A40                       CCAL0054
      CB240=CONS3*CA140                                                 CCAL0055
      CB340=8.*(A20*A21+2.*U*B40)/CONS1                                 CCAL0056
      C60=(CB140*ALPHA+4.*CB240*A20+16.*CB220*A40+2.*CB340*A21          CCAL0057
     1   +2.*CA120*B22+16.*CA220*B40+2.*CA340*B21+4.*CA320*B41          CCAL0058
     2   +2.*CC0*C42+4.*U*A21*C22)/(36.*CC1)                            CCAL0059
      CB141=4.*U  *C42+4.*A21*B22+2.*B41*ALPHA+16.*CONS3*A21*A40        CCAL0060
      CB341=8.*(A21**2+2.*A40*ALPHA+2.*U*B41)/CONS1                     CCAL0061
      C61=(4.*CA340*B22             +4.*CA321*B41+36.*CA201*B60         CCAL0062
     1   +16.*CA220*B41+CB141*ALPHA+4.*CB240*A21+16.*CB221*A40          CCAL0063
     2   +2.*CB341*A21+CA140*BETA+2.*CA121*B22+6.*CC0*C43+4.*U*C42*ALPHACCAL0064
     3   +12.*U*A21*C23)/(36.*CC1)                                      CCAL0065
      CB160=2.*U*C61+2.*A21*B41+16.*CONS3*A40**2                        CCAL0066
      CB360=8.*(2.*A21*A40+3.*U*B60)/CONS1                              CCAL0067
      C80=(CB160*ALPHA+16.*CB240*A40+2.*CB360*A21+2.*CA140*B22          CCAL0068
     1   +36.*CA220*B60+4.*CA340*B41+4.*U*A21*C42)/(64.*CC1)            CCAL0069
C                                                                       CCAL0070
 999  RETURN                                                            CCAL0071
      END                                                               CCAL0072
$IBFTC DCALC                                                            DCAL0000
      SUBROUTINE DCALC                                                  DCAL0001
      RETURN                                                            DCAL0002
      END                                                               DCAL0003
$IBFTC FCALC                                                            FCAL0000
CFCALC                                                                  FCAL0001
      SUBROUTINE FCALC                                                  FCAL0002
C                                                                       FCAL0003
        COMMON/ERR/IERR                                                 FCAL0004
      COMMON/NAMEA/VAR(6),IORDER,SOWCON,ISW1,ISW2                       FCAL0005
      COMMON/NAME1/THF,MARK,UUS,UDUDS,ZSINK                             FCAL0006
      COMMON/NAME2/VK,NFV,NVAR(6),X,U                                   FCAL0007
     1  ,ALPHA,BETA,GAMMA,DELTA,F(6),X2,X3,W(6)                         FCAL0008
     2  ,MF(6),D(6),AB(6,7),BLOCK(200),MAXIT,NUMIT                      FCAL0009
     3  ,X4,NSAV,XW,XW2,XW3,VKK,T2,XSTAR,Y                              FCAL0010
      DIMENSION A20(1),B21(1),C22(1),D23(1)                             FCAL0011
      COMMON/NAME3/                                                     FCAL0012
     1  A20,A21,A40,B21,B22,B40,B41,B42,B60,B61,B80,C22,C23,C40,C41,C42,FCAL0013
     2C43,C60,C61,C62,C63,C80,C81,C82,C100,C101,C120,D23,D24,D41,D42,D43FCAL0014
     3,D44,D60,D61,D62,D63,D64,D80,D81,D82,D83,D84,D100,D101,D102,D103  FCAL0015
     4 ,D120,D121,D122,D140,D141,D160,D40                               FCAL0016
     5 ,R,R2,R3,R4,R5,R6,R7,R8,R9,R10,R11,R12,R13,R14,R15               FCAL0017
     6 ,D1,D2,D3,D4,D5,D6,D7,D8,D9,D10,D11,D12,D13,D14,D15              FCAL0018
     7 ,UF,VF,DRDX1,DRDX2,DRDX3,DRDX4,DRDX5,DUDX1                       FCAL0019
     8 ,DUDX2,DUDX3,DUDX4,DVDX1,DVDX2,DVDX3,DVDX4                       FCAL0020
C                                                                       FCAL0021
      X=VAR(1)                                                          FCAL0022
      IF(SOWCON)33,31,33                                                FCAL0023
 31   IF(MARK)32,33,32                                                  FCAL0024
 32   TX=XW-X                                                           FCAL0025
      TX2=TX*TX                                                         FCAL0026
      TX3=TX*TX2                                                        FCAL0027
      VAR(2)=USTR+ASTR*TX+BSTR*TX2/2.0+GSTR*TX3/6.0                     FCAL0028
      VAR(3)=ASTR+BSTR*TX+GSTR*TX2/2.0                                  FCAL0029
 33   U=VAR(2)                                                          FCAL0030
      X2=X*X                                                            FCAL0031
      X3=X*X2                                                           FCAL0032
      X4=X*X3                                                           FCAL0033
C                                                                       FCAL0034
      BETA=0.0                                                          FCAL0035
      GAMMA=0.0                                                         FCAL0036
      DELTA=0.0                                                         FCAL0037
      GO TO (21,22,23,24),IORDER                                        FCAL0038
 24   DELTA=VAR(6)                                                      FCAL0039
 23   GAMMA=VAR(5)                                                      FCAL0040
 22   BETA=VAR(4)                                                       FCAL0041
 21   ALPHA=VAR(3)                                                      FCAL0042
      CALL ABCALC                                                       FCAL0043
      IF(IORDER-2) 4,4,2                                                FCAL0044
 2    CALL CCALC                                                        FCAL0045
      IF(IORDER-3) 4,4,3                                                FCAL0046
 3    CALL DCALC                                                        FCAL0047
 4    CALL PCALC                                                        FCAL0048
C                                                                       FCAL0049
      F(1)=UF                                                           FCAL0050
      F(2)=(U+UF)*DRDX1-VF                                              FCAL0051
      F(3)=DUDX1*DRDX1+(U+UF)*DRDX2-DVDX1                               FCAL0052
      F(4)=DUDX2*DRDX1+2.0*DUDX1*DRDX2+(U+UF)*DRDX3-DVDX2               FCAL0053
      F(5)=DUDX3*DRDX1+3.0*DUDX2*DRDX2+3.0*DUDX1*DRDX3                  FCAL0054
     1       +(U+UF)*DRDX4-DVDX3                                        FCAL0055
      F(6)=DUDX4*DRDX1+4.0*DUDX3*DRDX2+6.0*DUDX2*DRDX3                  FCAL0056
     1       +4.0*DUDX1*DRDX4+(U+UF)*DRDX5-DVDX4                        FCAL0057
      IF(SOWCON) 50,51,59                                               FCAL0058
 50   USTR=U+VAR(3)*X+VAR(4)*X2/2.0+VAR(5)*X3/6.0+VAR(6)*X4/24.0        FCAL0059
      ASTR=VAR(3)+VAR(4)*X+VAR(5)*X2/2.0+VAR(6)*X3/6.0                  FCAL0060
      GSTR=(4.0/ZSINK)*(USTR-UUS+(ASTR+UDUDS)*ZSINK/2.0)                FCAL0061
      BSTR=(UDUDS-ASTR-GSTR*(ZSINK**2)/2.0)/ZSINK                       FCAL0062
      GO TO 59                                                          FCAL0063
 51   IF(MARK) 52,59,52                                                 FCAL0064
 52   F(5)=0.0                                                          FCAL0065
      F(6)=0.0                                                          FCAL0066
      GO TO (53,53,54,59),IORDER                                        FCAL0067
 53   WRITE        (6,16)IORDER                                         FCAL0068
      IERR=1                                                            FCAL0069
 54   F(4)=0.0                                                          FCAL0070
C                                                                       FCAL0071
 59   DO 60 I=1,6                                                       FCAL0072
 60   F(I)=AMAX1(-1.E10,AMIN1(F(I),1.E10))                              FCAL0073
      IF(ISW1-1) 10,7,10                                                FCAL0074
 7    IF(NSAV-NUMIT)8,10,8                                              FCAL0075
 8    NSAV=NUMIT                                                        FCAL0076
      WRITE (6,11)NUMIT,(I,VAR(I),F(I),I=1,NFV)                         FCAL0077
      IF(ISW2-1) 10,9,10                                                FCAL0078
 9    WRITE (6,12)(A20(I),I=1,3),(B21(I),I=1,8),                        FCAL0079
     X(C22(I),I=1,16),(D23(I),I=1,26)                                   FCAL0080
 10   RETURN                                                            FCAL0081
C                                                                       FCAL0082
 11   FORMAT(17H0ITERATION NUMBER,I3,11X,1HI,6X,6HVAR(I),12X,           FCAL0083
     1   4HF(I),/(30X,I2,2E17.7))                                       FCAL0084
 12   FORMAT(36H0 1ST. ORDER COEFFICIENTS A20 TO A40,/3E15.7/,          FCAL0085
     1   36H  2ND. ORDER COEFFICIENTS B21 TO B80,/7E15.7/E15.7/,        FCAL0086
     2   37H  3RD. ORDER COEFFICIENTS C22 TO C120,/7E15.7/7E15.7/2E15.7/FCAL0087
     3   37H  4TH. ORDER COEFFICIENTS D23 TO D160,/(7E15.7))            FCAL0088
 16   FORMAT(1H0,28HERROR, CAN NOT FARE. IORDER=I2)                     FCAL0089
C                                                                       FCAL0090
      END                                                               FCAL0091
$IBFTC JAMES                                                            JAME0000
CJAMES                                                                  JAME0001
      SUBROUTINE JAMES (RR,XX,TT,VEK)                                   JAME0002
C                                                                       JAME0003
      EXTERNAL  FCALC                                                   JAME0004
C                                                                       JAME0005
        COMMON/ERR/IERR                                                 JAME0006
      COMMON/NAMEA/VAR(6),IORDER,SOWCON,ISW1,ISW2                       JAME0007
      COMMON/NAME1/THF,MARK,UUS,UDUDS,ZSINK                             JAME0008
      COMMON/NAME2/VK,NFV,NVAR(6),X,U                                   JAME0009
     1  ,ALPHA,BETA,GAMMA,DELTA,F(6),X2,X3,W(6)                         JAME0010
     2  ,MF(6),D(6),AB(6,7),BLOCK(200),MAXIT,NUMIT                      JAME0011
     3  ,X4,NSAV,XW,XW2,XW3,VKK,T2,XSTAR,Y                              JAME0012
      DIMENSION A20(1),B21(1),C22(1),D23(1)                             JAME0013
      COMMON/NAME3/                                                     JAME0014
     1  A20,A21,A40,B21,B22,B40,B41,B42,B60,B61,B80,C22,C23,C40,C41,C42,JAME0015
     2C43,C60,C61,C62,C63,C80,C81,C82,C100,C101,C120,D23,D24,D41,D42,D43JAME0016
     3,D44,D60,D61,D62,D63,D64,D80,D81,D82,D83,D84,D100,D101,D102,D103  JAME0017
     4 ,D120,D121,D122,D140,D141,D160,D40                               JAME0018
     5 ,R,R2,R3,R4,R5,R6,R7,R8,R9,R10,R11,R12,R13,R14,R15               JAME0019
     6 ,D1,D2,D3,D4,D5,D6,D7,D8,D9,D10,D11,D12,D13,D14,D15              JAME0020
     7 ,UF,VF,DRDX1,DRDX2,DRDX3,DRDX4,DRDX5,DUDX1                       JAME0021
     8 ,DUDX2,DUDX3,DUDX4,DVDX1,DVDX2,DVDX3,DVDX4                       JAME0022
C                                                                       JAME0023
      VK=VEK                                                            JAME0024
      VKK=(VK-1.0)/(VK+1.0)                                             JAME0025
      T2=1.0/TT**2                                                      JAME0026
      XW=XX                                                             JAME0027
      XW2=XW*XW                                                         JAME0028
      XW3=XW*XW2                                                        JAME0029
      DO 2 I=1,53                                                       JAME0030
 2    A20(I)=0.0                                                        JAME0031
C                                                                       JAME0032
      R=RR                                                              JAME0033
      R2=R*R                                                            JAME0034
      R3=R*R2                                                           JAME0035
      R4=R*R3                                                           JAME0036
      R5=R*R4                                                           JAME0037
      R6=R*R5                                                           JAME0038
      R7=R*R6                                                           JAME0039
      R8=R*R7                                                           JAME0040
      R9=R*R8                                                           JAME0041
      R10=R*R9                                                          JAME0042
      R11=R*R10                                                         JAME0043
      R12=R*R11                                                         JAME0044
      R13=R*R12                                                         JAME0045
      R14=R*R13                                                         JAME0046
      R15=R*R14                                                         JAME0047
C                                                                       JAME0048
      D1=2.0*R                                                          JAME0049
      D2=3.0*R2                                                         JAME0050
      D3=4.0*R3                                                         JAME0051
      D4=5.0*R4                                                         JAME0052
      D5=6.0*R5                                                         JAME0053
      D6=7.0*R6                                                         JAME0054
      D7=8.0*R7                                                         JAME0055
      D8=9.0*R8                                                         JAME0056
      D9=10.0*R9                                                        JAME0057
      D10=11.0*R10                                                      JAME0058
      D11=12.0*R11                                                      JAME0059
      D12=13.0*R12                                                      JAME0060
      D13=14.0*R13                                                      JAME0061
      D14=15.0*R14                                                      JAME0062
      D15=16.0*R15                                                      JAME0063
C                                                                       JAME0064
      MAXIT=15                                                          JAME0065
      NFV=IORDER+2                                                      JAME0066
      DO 1 I=1,NFV                                                      JAME0067
      MF(I)=0                                                           JAME0068
      NVAR(I)=0                                                         JAME0069
 1    W(I)=1.0                                                          JAME0070
C                                                                       JAME0071
      IF(SOWCON)20,23,21                                                JAME0072
 20   MARK=1                                                            JAME0073
      GO TO 22                                                          JAME0074
 21   MARK=0                                                            JAME0075
 22   MF(1)=1                                                           JAME0076
      NVAR(2)=1                                                         JAME0077
      GO TO 3                                                           JAME0078
 23   IF(XW+SIN(THF)/TT) 25,24,24                                       JAME0079
 24   MARK=0                                                            JAME0080
      GO TO 3                                                           JAME0081
 25   IF(MARK)26,3,26                                                   JAME0082
 26   NVAR(2)=1                                                         JAME0083
      NVAR(3)=1                                                         JAME0084
 3    NSAV=-1                                                           JAME0085
C                                                                       JAME0086
      CALL NEWT (NFV,F,MF,NFV,VAR,NVAR,FCALC,W                          JAME0087
     1  ,1.E-2,1.E-10,MAXIT,NUMIT,BLOCK,AB,D,0)                         JAME0088
C                                                                       JAME0089
      IF(IERR .NE. 0) GO TO 8                                           JAME0090
      IF (ISW1) 5,8,5                                                   JAME0091
 8    WRITE (6,19)NUMIT,(I,VAR(I),F(I),I=1,NFV)                         JAME0092
 5    IF(NUMIT-MAXIT) 11,9,9                                            JAME0093
 9    WRITE (6,18)MAXIT                                                 JAME0094
      IERR=1                                                            JAME0095
 11   RETURN                                                            JAME0096
C                                                                       JAME0097
 18   FORMAT(22H0NO CONVERGENCE WITHINI3,1X,10HITERATIONS1X,///)        JAME0098
 19   FORMAT(17H0ITERATION NUMBER,I3,11X,1HI,6X,6HVAR(I),12X,           JAME0099
     1   4HF(I),/(30X,I2,2E17.7))                                       JAME0100
      END                                                               JAME0101
$IBFTC LEGS                                                             LEGS0000
      SUBROUTINE LEGS(AB,DPAR,BNDS,M,N1,MD,ND,ATA,R,Z,                  LEGS0001
     XSUS,SUSP,I1,I2,I3,I4,I5,IDLE)                                     LEGS0002
C REMOVING THE C FROM THE NEXT TWO CARDS MAKES LEGS DOUBLE PRECISION    LEGS0003
C     DOUBLE PRECISION AB,DPAR,BNDS,ATA,R,Z,SUS,SUSP,D1,D2,D            LEGS0004
C    1,E1,E2,H,HAL,RHI,RHV,S1,S,SS,X1,X2,X3,X,Y1,Y2,Y3,Y,X1AL,TP        LEGS0005
      DIMENSION AB(2),DPAR(2),BNDS(2),ATA(2),R(2),IDLE(2)               LEGS0006
                                                                        LEGS0007
C  COUNT THE NUMBER OF DELETIONS                                        LEGS0008
                                                                        LEGS0009
      NIDLE=0                                                           LEGS0010
      J=1                                                               LEGS0011
    1 IF(IDLE(J))3,3,2                                                  LEGS0012
    2 NIDLE=NIDLE+1                                                     LEGS0013
      J=J+2                                                             LEGS0014
      IF(J-2*M)1,1,3                                                    LEGS0015
                                                                        LEGS0016
C  SET UP COUNTERS                                                      LEGS0017
                                                                        LEGS0018
    3 N=N1                                                              LEGS0019
      N2=N+1                                                            LEGS0020
      IQ=(N*N2)/2                                                       LEGS0021
      IP=(N2*(N2+1))/2                                                  LEGS0022
                                                                        LEGS0023
                                                                        LEGS0024
C  BRANCH ON I3                                                         LEGS0025
                                                                        LEGS0026
C     I3 LESS THAN ZERO -- COMPUTE ATA AND ADD IT TO CONTENTS OF ATA    LEGS0027
C     I3 EQUAL TO ZERO  -- COMPUTE ATA AND PUT IN CLEARED ATA           LEGS0028
C     I3 GREATER THAN ZERO -- ATA CALCULATED BEFORE ENTRY               LEGS0029
                                                                        LEGS0030
                                                                        LEGS0031
      IF(I3)6,4,14                                                      LEGS0032
    4 DO 5 I=1,IP                                                       LEGS0033
    5 ATA(I)=0.                                                         LEGS0034
    6 DO 13 K=1,M                                                       LEGS0035
      IF(NIDLE)10,10,7                                                  LEGS0036
                                                                        LEGS0037
C  ROW DELETOR SECTION                                                  LEGS0038
                                                                        LEGS0039
    7 DO 9 J=1,NIDLE                                                    LEGS0040
      IS=2*J                                                            LEGS0041
      IF(K-IDLE(IS-1))9,13,8                                            LEGS0042
    8 IF(K-IDLE(IS))13,13,9                                             LEGS0043
    9 CONTINUE                                                          LEGS0044
                                                                        LEGS0045
C  COMPUTE ATA                                                          LEGS0046
                                                                        LEGS0047
   10 L=1                                                               LEGS0048
      LI=K                                                              LEGS0049
      DO 12 I=1,N2                                                      LEGS0050
      LJ=LI                                                             LEGS0051
      DO 11 J=I,N2                                                      LEGS0052
      ATA(L)=ATA(L)+AB(LI)*AB(LJ)                                       LEGS0053
      L=L+1                                                             LEGS0054
   11 LJ=LJ+MD                                                          LEGS0055
   12 LI=LI+MD                                                          LEGS0056
   13 CONTINUE                                                          LEGS0057
                                                                        LEGS0058
                                                                        LEGS0059
C BRANCH ON I5                                                          LEGS0060
                                                                        LEGS0061
C     I5 EQUAL TO ZERO     -- CONTINUE                                  LEGS0062
C     I5 NOT EQUAL TO ZERO -- EXIT                                      LEGS0063
                                                                        LEGS0064
                                                                        LEGS0065
      IF(I5)82,14,82                                                    LEGS0066
                                                                        LEGS0067
                                                                        LEGS0068
C  BRANCH ON I2                                                         LEGS0069
                                                                        LEGS0070
C     I2 EQUAL TO ZERO     -- SOLVE FOR X                               LEGS0071
C     I2 NOT EQUAL TO ZERO -- DO NOT SOLVE FOR X                        LEGS0072
                                                                        LEGS0073
                                                                        LEGS0074
   14 IF(I2)15,17,15                                                    LEGS0075
   15 DO 16 I=1,N                                                       LEGS0076
   16 DPAR(I)=BNDS(I)                                                   LEGS0077
      X=1.                                                              LEGS0078
      GO TO 71                                                          LEGS0079
                                                                        LEGS0080
                                                                        LEGS0081
C  BRANCH ON I1                                                         LEGS0082
                                                                        LEGS0083
C     I1 EQUAL TO ZERO     -- NO BOUNDS IN SOLVING FOR X                LEGS0084
C     I1 NOT EQUAL TO ZERO -- BOUNDS USED IN SOLVING FOR X              LEGS0085
                                                                        LEGS0086
   17 IF(I1)20,18,20                                                    LEGS0087
   18 DO 19 I=1,N                                                       LEGS0088
   19 BNDS(I)=-1.                                                       LEGS0089
                                                                        LEGS0090
C  INITIALIZING FOR ITERATION                                           LEGS0091
                                                                        LEGS0092
   20 RHV=.0001                                                         LEGS0093
      RHI=.2                                                            LEGS0094
      E1=.1                                                             LEGS0095
      E2=.2                                                             LEGS0096
      IMAX=20                                                           LEGS0097
      ICT=0                                                             LEGS0098
      X=0.                                                              LEGS0099
      IF=1                                                              LEGS0100
                                                                        LEGS0101
C  FIND UNBOUNDED SOLUTION                                              LEGS0102
                                                                        LEGS0103
      GO TO 47                                                          LEGS0104
                                                                        LEGS0105
C  IF CONSTRAINTS SATISFIED,EXIT ITERATION                              LEGS0106
                                                                        LEGS0107
   21 IF(Y-E1)71,71,22                                                  LEGS0108
                                                                        LEGS0109
C  GET FIRST ESTIMATE FOR ITERATION                                     LEGS0110
                                                                        LEGS0111
   22 H=0.                                                              LEGS0112
      KT=1                                                              LEGS0113
      KL=N2                                                             LEGS0114
      IF=2                                                              LEGS0115
      K=KT+KL-1                                                         LEGS0116
      HAL=0.0                                                           LEGS0117
      DO 28 I=1,N                                                       LEGS0118
      IF (BNDS(I)) 27,27,23                                             LEGS0119
 23   X1=DABS(ATA(K))                                                   LEGS0120
      X1=BNDS(I)**2*(X1/BNDS(I)-ATA(KT))                                LEGS0121
      X1AL=ATA(KT)*BNDS(I)**2*RHV                                       LEGS0122
      IF(X1AL-HAL)25,25,24                                              LEGS0123
   24 HAL=X1AL                                                          LEGS0124
   25 IF(X1-H)27,27,26                                                  LEGS0125
   26 H=X1                                                              LEGS0126
   27 KT=KT+KL                                                          LEGS0127
      KL=KL-1                                                           LEGS0128
   28 K=K+KL                                                            LEGS0129
      IF(H)29,29,30                                                     LEGS0130
   29 H=HAL                                                             LEGS0131
                                                                        LEGS0132
C  SET X1=FIRST ESTIMATE                                                LEGS0133
                                                                        LEGS0134
   30 X1=X                                                              LEGS0135
      Y1=Y                                                              LEGS0136
      X=X+H                                                             LEGS0137
      GO TO 47                                                          LEGS0138
   31 IF(Y-E1)33,33,32                                                  LEGS0139
                                                                        LEGS0140
C  X IS TOO SMALL INCREASE AND TRY AGAIN                                LEGS0141
                                                                        LEGS0142
   32 H=10.*H                                                           LEGS0143
      GO TO 30                                                          LEGS0144
                                                                        LEGS0145
C  SOLUTION IS BETWEEN X1 AND X.                                        LEGS0146
C  IF X IS GOOD ENOUGH, EXIT ITERATION                                  LEGS0147
                                                                        LEGS0148
   33 IF(Y+E2)34,71,71                                                  LEGS0149
   34 X2=X                                                              LEGS0150
      Y2=Y                                                              LEGS0151
      IF=3                                                              LEGS0152
                                                                        LEGS0153
C  SOLUTION IS BETWEEN X1 AND X2. TRY A POINT IN BETWEEN.               LEGS0154
                                                                        LEGS0155
   35 X=(1.-RHI)*X1+RHI*X2                                              LEGS0156
      IF=3                                                              LEGS0157
      GO TO 47                                                          LEGS0158
   36 IF(Y-E1)37,37,39                                                  LEGS0159
                                                                        LEGS0160
C  IF X IS GOOD ENOUGH, EXIT ITERATION                                  LEGS0161
                                                                        LEGS0162
   37 IF(Y+E1)38,71,71                                                  LEGS0163
                                                                        LEGS0164
C  SOLUTION IS BETWEEN X1 AND X.                                        LEGS0165
                                                                        LEGS0166
   38 X3=X2                                                             LEGS0167
      Y3=Y2                                                             LEGS0168
      X2=X                                                              LEGS0169
      Y2=Y                                                              LEGS0170
                                                                        LEGS0171
C  SOLUTION IS BETWEEN X1 AND X2. X3 BIGGER THAN X2                     LEGS0172
                                                                        LEGS0173
      GO TO 40                                                          LEGS0174
                                                                        LEGS0175
C  SOLUTION IS BETWEEN X AND X2                                         LEGS0176
                                                                        LEGS0177
   39 X3=X1                                                             LEGS0178
      Y3=Y1                                                             LEGS0179
      X1=X                                                              LEGS0180
      Y1=Y                                                              LEGS0181
                                                                        LEGS0182
C  SOLUTION IS BETWEEN X1 AND X2. X3 LESS THAN X1                       LEGS0183
                                                                        LEGS0184
   40 D1=Y2-Y1                                                          LEGS0185
                                                                        LEGS0186
C  DO INVERSE QUADRATIC INTERPOLATION TO GET NEXT POINT                 LEGS0187
                                                                        LEGS0188
      D2=Y3-Y2                                                          LEGS0189
      D=Y3-Y1                                                           LEGS0190
      X=0.                                                              LEGS0191
      TP=D1*D                                                           LEGS0192
      IF(TP)401,402,401                                                 LEGS0193
 401  X=X+(X1*Y2*Y3)/TP                                                 LEGS0194
 402  TP=D1*D2                                                          LEGS0195
      IF(TP)403,404,403                                                 LEGS0196
 403  X=X-(X2*Y1*Y3)/TP                                                 LEGS0197
 404  TP=D*D2                                                           LEGS0198
      IF(TP)405,406,405                                                 LEGS0199
 405  X=X+(X3*Y1*Y2)/TP                                                 LEGS0200
 406  CONTINUE                                                          LEGS0201
                                                                        LEGS0202
C  IF INTERPOLATED POINT IS OUTSIDE RANGE (X1,X2),                      LEGS0203
C  TRY A NEW POINT BETWEEN X1 AND X2                                    LEGS0204
                                                                        LEGS0205
      IF(X-X2)41,35,35                                                  LEGS0206
   41 IF(X-X1)35,35,42                                                  LEGS0207
                                                                        LEGS0208
C  INTERPOLATED POINT IS WITHIN RANGE (X1,X2)                           LEGS0209
                                                                        LEGS0210
   42 IF=4                                                              LEGS0211
      GO TO 47                                                          LEGS0212
                                                                        LEGS0213
C  IF X IS GOOD ENOUGH, EXIT ITERATION                                  LEGS0214
                                                                        LEGS0215
   43 IF(Y-E1)44,44,46                                                  LEGS0216
   44 IF(Y+E2)45,71,71                                                  LEGS0217
   45 X2=X                                                              LEGS0218
                                                                        LEGS0219
C  SOLUTION IS BETWEEN X1 AND NEW X2                                    LEGS0220
                                                                        LEGS0221
      Y2=Y                                                              LEGS0222
      GO TO 35                                                          LEGS0223
   46 X1=X                                                              LEGS0224
      Y1=Y                                                              LEGS0225
                                                                        LEGS0226
C  SOLUTION IS BETWEEN NEW X1 AND X2                                    LEGS0227
                                                                        LEGS0228
      GO TO 35                                                          LEGS0229
                                                                        LEGS0230
C  SOLVING THE LINEAR SYSTEM                                            LEGS0231
C  DECOMPOSE X SUBMATRIX                                                LEGS0232
                                                                        LEGS0233
   47 R(1)=-1.                                                          LEGS0234
      R(IP+1)=0.                                                        LEGS0235
      L=1                                                               LEGS0236
      IF(BNDS(1))49,50,48                                               LEGS0237
   48 R(IP+1)=X/BNDS(1)**2                                              LEGS0238
                                                                        LEGS0239
C  COMPUTE FIRST ELEMENT OF D.                                          LEGS0240
                                                                        LEGS0241
   49 R(IP+1)=ATA(1)+R(IP+1)                                            LEGS0242
   50 DO 62 K=1,N                                                       LEGS0243
                                                                        LEGS0244
C  DECOMPOSE (K+1)X(K+1) SUBMATRIX                                      LEGS0245
                                                                        LEGS0246
      IF(K-N)51,53,53                                                   LEGS0247
   51 IF(BNDS(K+1))53,52,53                                             LEGS0248
   52 K1=IP+K+1                                                         LEGS0249
      R(K1)=0.                                                          LEGS0250
      GO TO 61                                                          LEGS0251
   53 KL=1                                                              LEGS0252
      DO 55 J=1,K                                                       LEGS0253
      KM=K+1                                                            LEGS0254
      KT=N                                                              LEGS0255
      S=0.                                                              LEGS0256
      DO 54 I=1,J                                                       LEGS0257
      S=S+R(KL)*ATA(KM)                                                 LEGS0258
      KM=KM+KT                                                          LEGS0259
      KT=KT-1                                                           LEGS0260
   54 KL=KL+1                                                           LEGS0261
      K1=L+J                                                            LEGS0262
      K2=IP+J                                                           LEGS0263
      R(K1)=0.                                                          LEGS0264
      IF(R(K2))541,55,541                                               LEGS0265
 541  R(K1)=S/R(K2)                                                     LEGS0266
 55   CONTINUE                                                          LEGS0267
      KL=0                                                              LEGS0268
      DO 57 J=1,K                                                       LEGS0269
      KL1=KL+J                                                          LEGS0270
      S=0.                                                              LEGS0271
      DO 56 I=J,K                                                       LEGS0272
      K1=L+I                                                            LEGS0273
      S=S+R(KL1)*R(K1)                                                  LEGS0274
   56 KL1=KL1+I                                                         LEGS0275
      K1=L+J                                                            LEGS0276
      KL=KL+J                                                           LEGS0277
   57 R(K1)=S                                                           LEGS0278
      KL=K+1                                                            LEGS0279
      KT=N                                                              LEGS0280
      S=0.                                                              LEGS0281
      DO 58 J=1,K                                                       LEGS0282
      K1=L+J                                                            LEGS0283
      S=S+R(K1)*ATA(KL)                                                 LEGS0284
      KL=KL+KT                                                          LEGS0285
   58 KT=KT-1                                                           LEGS0286
      K1=IP+K+1                                                         LEGS0287
      R(K1)=ATA(KL)-S                                                   LEGS0288
      IF(K-N)59,61,61                                                   LEGS0289
   59 IF(BNDS(K+1))61,61,60                                             LEGS0290
   60 R(K1)=R(K1)+X/(BNDS(K+1))**2                                      LEGS0291
   61 L=L+K+1                                                           LEGS0292
      R(L)=-1.                                                          LEGS0293
   62 CONTINUE                                                          LEGS0294
      GO TO(63,63,63,63,79),IF                                          LEGS0295
   63 DO 64 I=1,N                                                       LEGS0296
      K1=IQ+I                                                           LEGS0297
   64 DPAR(I)=R(K1)                                                     LEGS0298
      ICT=ICT+1                                                         LEGS0299
      IF(ICT-IMAX)65,71,71                                              LEGS0300
   65 S=0.                                                              LEGS0301
      DO 67 I=1,N                                                       LEGS0302
      IF(BNDS(I))67,67,66                                               LEGS0303
   66 S=S+(DPAR(I)/BNDS(I))**2                                          LEGS0304
   67 CONTINUE                                                          LEGS0305
   68 Y=S-1.                                                            LEGS0306
   69 Z=X                                                               LEGS0307
   70 GO TO(21,31,36,43,79),IF                                          LEGS0308
                                                                        LEGS0309
C  FIND LENGTH OF RESIDUAL VECTOR                                       LEGS0310
                                                                        LEGS0311
   71 SS=0.                                                             LEGS0312
      S1=0.                                                             LEGS0313
      K1=N2                                                             LEGS0314
      K2=N                                                              LEGS0315
      DO 76 J=1,N                                                       LEGS0316
      L=J                                                               LEGS0317
      S=0.                                                              LEGS0318
      KK=J-1                                                            LEGS0319
      IF(KK)74,74,72                                                    LEGS0320
   72 DO 73 I=1,KK                                                      LEGS0321
      S=S+ATA(L)*DPAR(I)                                                LEGS0322
   73 L=L+N2-I                                                          LEGS0323
   74 DO 75 I=J,N                                                       LEGS0324
      S=S+ATA(L)*DPAR(I)                                                LEGS0325
   75 L=L+1                                                             LEGS0326
      SS=SS+S*DPAR(J)                                                   LEGS0327
      S1=S1+ATA(K1)*DPAR(J)                                             LEGS0328
      K1=K1+K2                                                          LEGS0329
   76 K2=K2-1                                                           LEGS0330
      SUS=ATA(K1)                                                       LEGS0331
      SUSP=SS-2.*S1+SUS                                                 LEGS0332
                                                                        LEGS0333
                                                                        LEGS0334
C  BRANCH ON I4                                                         LEGS0335
                                                                        LEGS0336
C     I4 EQUAL TO ZERO     -- INVERT NORMAL MATRIX                      LEGS0337
C     I4 NOT EQUAL TO ZERO -- DO NOT INVERT NORMAL MATRIX,EXIT          LEGS0338
                                                                        LEGS0339
      IF(I4)82,77,82                                                    LEGS0340
   77 IF(X)78,79,78                                                     LEGS0341
   78 X=0.                                                              LEGS0342
      IF=5                                                              LEGS0343
      GO TO 47                                                          LEGS0344
                                                                        LEGS0345
C  FIND THE INVERSE OF (A TRANSPOSE)*A                                  LEGS0346
                                                                        LEGS0347
   79 IRSS=0                                                            LEGS0348
      DO 81 K=1,N                                                       LEGS0349
      IRSS=IRSS+K                                                       LEGS0350
      IRS=IRSS-K                                                        LEGS0351
      ITS=IRSS                                                          LEGS0352
      DO 81 J=K,N                                                       LEGS0353
      IRS=IRS+J                                                         LEGS0354
      IT=ITS                                                            LEGS0355
      IR=IRS                                                            LEGS0356
      S=0.                                                              LEGS0357
      DO 80 I=J,N                                                       LEGS0358
      IN=IP+I                                                           LEGS0359
      IF(R(IN))791,792,791                                              LEGS0360
 791  S=S+(R(IT)*R(IR))/R(IN)                                           LEGS0361
 792  IT=IT+1                                                           LEGS0362
   80 IR=IR+I                                                           LEGS0363
      R(ITS)=S                                                          LEGS0364
   81 ITS=ITS+J                                                         LEGS0365
   82 RETURN                                                            LEGS0366
      END                                                               LEGS0367
$IBFTC NEWT                                                             NEWT0000
      SUBROUTINE NEWT(MM,F,MF,NN,VAR,NVAR,FCALC,WF,                     NEWT0001
     1 PERTV,EPS1,MAXIT,NUMIT,BIN,AB,B,KB)                              NEWT0002
C                                                                       NEWT0003
C      KB=0, UNBOUNDED                                                  NEWT0004
C      KB=1, BOUNDED                                                    NEWT0005
C                                                                       NEWT0006
      LOGICAL TF                                                        NEWT0007
      DIMENSION F(MM),MF(MM),VAR(NN),NVAR(NN),WF(MM),                   NEWT0008
     1 AB(MM,1),BIN(1),B(NN)                                            NEWT0009
C                                                                       NEWT0010
C     SET-UP CALCULATIONS                                               NEWT0011
C                                                                       NEWT0012
      M=MM                                                              NEWT0013
      N=NN                                                              NEWT0014
      NUMIT=0                                                           NEWT0015
      TF=.FALSE.                                                        NEWT0016
      NVPL1=N+1                                                         NEWT0017
      NR=3*N+((N+1)*(N+2))/2+1                                          NEWT0018
      NUMF=0                                                            NEWT0019
      DO 10 I=1,M                                                       NEWT0020
      IF(MF(I).NE.0) GO TO 10                                           NEWT0021
      NUMF=NUMF+1                                                       NEWT0022
 10   CONTINUE                                                          NEWT0023
      NUMVAR=0                                                          NEWT0024
      N2=2*N                                                            NEWT0025
      J=N2                                                              NEWT0026
      DO 13 I=1,N                                                       NEWT0027
      IF(NVAR(I).NE.0) GO TO 13                                         NEWT0028
      NUMVAR=NUMVAR+1                                                   NEWT0029
      BIN(I)=VAR(I)                                                     NEWT0030
      J=J+1                                                             NEWT0031
      BIN(J)=B(I)                                                       NEWT0032
 13   CONTINUE                                                          NEWT0033
      FSUML=1.E38                                                       NEWT0034
C                                                                       NEWT0035
C     ITERATION LOOP                                                    NEWT0036
C                                                                       NEWT0037
 14   DO 29 L=1,MAXIT                                                   NEWT0038
      CALL FCALC                                                        NEWT0039
      FSUM=0.                                                           NEWT0040
      J=0                                                               NEWT0041
      DO 16 I=1,M                                                       NEWT0042
      IF(MF(I).NE.0) GO TO 16                                           NEWT0043
      J=J+1                                                             NEWT0044
      AB(J,NUMVAR+1)=-F(I)                                              NEWT0045
      FSUM=FSUM+(WF(I)*F(I))**2                                         NEWT0046
 16   CONTINUE                                                          NEWT0047
      IF(FSUM.LE.EPS1)RETURN                                            NEWT0048
      IF(TF) RETURN                                                     NEWT0049
      IF(FSUML.GT.FSUM) GO TO 22                                        NEWT0050
C                                                                       NEWT0051
C      HALVE BACK                                                       NEWT0052
      DO 21 I=1,N                                                       NEWT0053
      IF(NVAR(I).NE.0) GO TO 21                                         NEWT0054
      VAR(I)=(BIN(I)+VAR(I))/2.                                         NEWT0055
 21   CONTINUE                                                          NEWT0056
      GO TO 29                                                          NEWT0057
C                                                                       NEWT0058
 22   FSUML=FSUM                                                        NEWT0059
C                                                                       NEWT0060
C     CALCULATE PARTIAL DERIVATIVES                                     NEWT0061
C                                                                       NEWT0062
      J1=0                                                              NEWT0063
      DO 26 J=1,N                                                       NEWT0064
      IF(NVAR(J).NE.0) GO TO 26                                         NEWT0065
      J1=J1+1                                                           NEWT0066
      PERTV1=AMAX1(ABS(VAR(J)*PERTV),1.E-3*PERTV)                       NEWT0067
      VAR(J)=VAR(J)+PERTV1                                              NEWT0068
      CALL FCALC                                                        NEWT0069
      VAR(J)=VAR(J)-PERTV1                                              NEWT0070
      II=0                                                              NEWT0071
      DO 25 I=1,M                                                       NEWT0072
      IF(MF(I).NE.0) GO TO 25                                           NEWT0073
      II=II+1                                                           NEWT0074
      AB(II,J1)=(F(I)+AB(II,NUMVAR+1))/PERTV1                           NEWT0075
 25   CONTINUE                                                          NEWT0076
 26   CONTINUE                                                          NEWT0077
C                                                                       NEWT0078
C      CALCULATE NEXT ESTIMATE OF SOLUTION VECTOR                       NEWT0079
      CALL LEGS(AB,BIN(NVPL1),BIN(N2+1),NUMF,NUMVAR,                    NEWT0080
     1 M,NVPL1,BIN(3*N+1),BIN(NR),Z,SUS,SUSP,KB,0,0,0,0,0)              NEWT0081
      I=N                                                               NEWT0082
      DO 28 J=1,N                                                       NEWT0083
      IF(NVAR(J).NE.0) GO TO 28                                         NEWT0084
      I=I+1                                                             NEWT0085
      BIN(J)=VAR(J)                                                     NEWT0086
      VAR(J)=VAR(J)+BIN(I)                                              NEWT0087
 28   CONTINUE                                                          NEWT0088
 29   NUMIT=L                                                           NEWT0089
      TF=.TRUE.                                                         NEWT0090
      GO TO 14                                                          NEWT0091
C                                                                       NEWT0092
      END                                                               NEWT0093
$IBFTC ONED                                                             ONED0000
      SUBROUTINE ONED                                                   ONED0001
C                                                                       ONED0002
      DIMENSION  WPWG(10),PS(50,50),PW(150),R(10),XXXX(4881),FA(100),   ONED0003
     XGF(10),WDP(10),CK(10),CL(10),CKP(10),WPWT(10),SAV(6,10),          ONED0004
     XSAVE(53,10),R234(10),XXX(20),                                     ONED0005
     XUP(10),TP(10),DUP(10),DTP(10)                                     ONED0006
      COMMON CAPN  ,CPL   ,CPS   ,CTOP  ,DL    ,DR    ,EK    ,EKW   ,   ONED0007
     XEW    ,WPWG  ,GAMMA ,GMG0  ,HPL   ,HPS   ,IMAX  ,N1    ,N2    ,   ONED0008
     XPS    ,PG0   ,PW    ,R     ,RCAP  ,RHOP0 ,SMP   ,TG0   ,TPM   ,   ONED0009
     XSAV   ,SAVE  ,XXXX  ,FA    ,VGR   ,CPG   ,NILP  ,THI   ,GF    ,   ONED0010
     XWDG   ,IFG   ,RHOG0 ,THRUST,TI    ,WDP   ,TMD   ,IPPX  ,BC    ,   ONED0011
     XTHII  ,DZMIN ,ITOT  ,ZAX   ,RRT   ,RT    ,VK    ,T     ,CG    ,   ONED0012
     XA0C   ,DZI   ,ZI    ,ZJ    ,NZONES,CK    ,CL    ,CKP   ,R1    ,   ONED0013
     XZ1    ,UG    ,VG    ,TT0   ,DG    ,UGR   ,UGZ   ,VRG   ,VGZ   ,   ONED0014
     XZW    ,R234  ,WPWGT ,WPWT  ,XXX   ,USG                            ONED0015
      COMMON/NAME1/THF,MARK,UUS,UDUDS,ZSINK                             ONED0016
C                                                                       ONED0017
C         SET OF INITIAL CONDITIONS                                     ONED0018
C                                                                       ONED0019
      B=(1.+WPWGT)/(1.+WPWGT*CPL/CPG)                                   ONED0020
      G1=GAMMA-1.                                                       ONED0021
      C=1.+WPWGT*(1.+G1*CPL/CPG*B)                                      ONED0022
      UGMAX=SQRT (2.*CPG*TG0/B)                                         ONED0023
      G2=G1*B/C                                                         ONED0024
      GBAR=G2+1.                                                        ONED0025
      Q1=199.+GBAR                                                      ONED0026
      UG=SQRT (G2/Q1)*UGMAX                                             ONED0027
      Q2=200./Q1                                                        ONED0028
      TPI0=Q2*TG0                                                       ONED0029
      DO 3 I=1,ITOT                                                     ONED0030
      UP(I)=UG                                                          ONED0031
 3    TP(I)=TPI0                                                        ONED0032
      A=WDG/Q2**(1./G2)/UG/RHOG0                                        ONED0033
      C1THI=1.-COS (THI)                                                ONED0034
      AC=6.2831853*RT*RT*C1THI                                          ONED0035
      Z=-SQRT (A/AC)                                                    ONED0036
      ZMAX=-(1.+RRT*C1THI)/SIN (THI)                                    ONED0037
      DZ=DZMIN                                                          ONED0038
      WRITE (6,90)Z,ZMAX,DZ,UG,TPI0,A                                   ONED0039
 90   FORMAT(21H0ONE DIMENSIONAL FLOW/4H Z0=F9.4,                       ONED0040
     X 4X,4HMAX=F9.5,4X,10HSTEP SIZE=F9.6,4X,3HU0=F9.3,                 ONED0041
     X 4X,5HTPI0=F9.3,4X,3HA0=F9.4/)                                    ONED0042
C                                                                       ONED0043
C         CONSTANTS FOR INTEGRATION                                     ONED0044
C                                                                       ONED0045
      CUT=-A0C/1.5*CTOP/CPL                                             ONED0046
      CU1=GAMMA*RCAP                                                    ONED0047
      CU2=RCAP+RCAP                                                     ONED0048
      CU4=RCAP/CPG                                                      ONED0049
      CU3=CU4*CPL                                                       ONED0050
C                                                                       ONED0051
C         INTEGRATION WITH PRINT SELECTION                              ONED0052
C                                                                       ONED0053
      TMDR=TMD+Z                                                        ONED0054
      LT=0                                                              ONED0055
      ASSIGN 13 TO LTT                                                  ONED0056
      DO 15 JJ=1,10000                                                  ONED0057
      IF (Z+DZ-ZMAX) 5,4,4                                              ONED0058
 4    LT=1                                                              ONED0059
      ASSIGN 14 TO LTT                                                  ONED0060
      DZ=ZMAX-Z                                                         ONED0061
 5    A=AC*Z**2                                                         ONED0062
      RHOG=WDG/UG/A                                                     ONED0063
      T1=0.                                                             ONED0064
      DO 6 I=1,ITOT                                                     ONED0065
 6    T1=T1+WPWG(I)*(CPL*(TP(I)-TG0)+UP(I)**2/2.)                       ONED0066
      TG=TG0-(UG**2/2.+T1)/CPG                                          ONED0067
      TT0N=(TG/TG0)**CAPN                                               ONED0068
      Q4=SQRT (RCAP*TG)                                                 ONED0069
      GMG=GMG0*TT0N                                                     ONED0070
      DU1=0.                                                            ONED0071
      DU2=0.                                                            ONED0072
      DO 10 I=1,ITOT                                                    ONED0073
      Q3=GMG/RHOG/R(I)                                                  ONED0074
      CBAR=Q3/Q4                                                        ONED0075
      RE=2.*(UG-UP(I))/Q3                                               ONED0076
      CALL FIND(IFG,FA,RE,FP)                                           ONED0077
      F=FP*(1.+2.52*CBAR)/(1.+3.78*CBAR)                                ONED0078
      G=FP/(1.+CG*CBAR*FP)                                              ONED0079
      DUT=TT0N/R(I)**2/UP(I)                                            ONED0080
      DUP(I)=A0C*DUT*(UG-UP(I))*F                                       ONED0081
      DTP(I)=CUT*DUT*(TP(I)-TG)*G                                       ONED0082
      DU1=WPWG(I)*DTP(I)+DU1                                            ONED0083
      DU2=WPWG(I)*(UG-CU4*UP(I))*DUP(I)+DU2                             ONED0084
      UP(I)=UP(I)+DUP(I)*DZ                                             ONED0085
 10   TP(I)=TP(I)+DTP(I)*DZ                                             ONED0086
      DUG=UG*GAMMA/(UG**2-CU1*TG)*(CU2*TG/Z+CU3*DU1-DU2)                ONED0087
      UG=UG+DUG*DZ                                                      ONED0088
      Z=Z+DZ                                                            ONED0089
      GO TO LTT, (13,14)                                                ONED0090
 13   IF (Z-TMDR) 15,15,14                                              ONED0091
 14   TMDR=TMDR+TMD                                                     ONED0092
      AAT=2.*C1THI*Z**2                                                 ONED0093
      DO 12 I=1,ITOT                                                    ONED0094
      CK(I)=UP(I)/UG                                                    ONED0095
 12   CL(I)=(TG0-TP(I))/(TG0-TG)                                        ONED0096
      WRITE (6,91)Z,AAT,UG,ITOT,                                        ONED0097
     X(UP(I),TP(I),CK(I),CL(I),I=1,ITOT)                                ONED0098
 91   FORMAT (3H0Z=F9.5,5X5HA/AT=F9.5,8X3HUG=F9.3,5X5HUP(I)5X5HTP(I)    ONED0099
     X 8X,4HK(I),6X,4HL(I),5X,6HI=1,..I3/(2F12.3,2F12.5,2F12.3,2F12.5)) ONED0100
      IF(LT) 20,15,20                                                   ONED0101
 15   CONTINUE                                                          ONED0102
 20   UUS=(UG-DUG*DZ)/USG                                               ONED0103
      UDUDS=DUG/USG                                                     ONED0104
      STHI=SIN (THI)                                                    ONED0105
      CTHI=COS (THI)                                                    ONED0106
      RBAR=(1.+RRT*(1.-CTHI))/STHI                                      ONED0107
      ZSINK=Z-DZ+RBAR*CTHI-RRT*STHI                                     ONED0108
      RETURN                                                            ONED0109
       END                                                              ONED0110
$IBFTC PARTIL                                                           PART0000
      SUBROUTINE PARTIL                                                 PART0001
C                                                                       PART0002
      DIMENSION  WPWG(10),PS(50,50),PW(150),R(10),XXXX(4881),FA(100),   PART0003
     XGF(10),WDP(10),CK(10),CL(10),CKP(10),WPWT(10),SAV(6,10),          PART0004
     XSAVE(53,10),R234(10),CKK(10),CLL(10),ZW1(10),SAUR(6)              PART0005
     X,R2(10,4),                                                        PART0006
     XCKU0(10),CLT0(10),UP(10),VP(10),TP(10),PSI(10,4),                 PART0007
     XCKU(10,4),CKV(10,4),CLT(10,4),                                    PART0008
     XT1(10),T2(10),T3(10),T4(10),T5(10),T6(10),T7(10)                  PART0009
      COMMON CAPN  ,CPL   ,CPS   ,CTOP  ,DL    ,DR    ,EK    ,EKW   ,   PART0010
     XEW    ,WPWG  ,GAMMA ,GMG0  ,HPL   ,HPS   ,IMAX  ,N1    ,N2    ,   PART0011
     XPS    ,PG0   ,PW    ,R     ,RCAP  ,RHOP0 ,SMP   ,TG0   ,TPM   ,   PART0012
     XSAV   ,SAVE  ,XXXX  ,FA    ,VGR   ,CPG   ,NILP  ,THI   ,GF    ,   PART0013
     XWDG   ,IFG   ,RHOG0 ,THRUST,TI    ,WDP   ,TMD   ,IPPX  ,BC    ,   PART0014
     XNIL   ,DZMIN ,ITOT  ,ZAX   ,RRT   ,RT    ,VK    ,T     ,CG    ,   PART0015
     XA0C   ,DZI   ,ZI    ,ZJ    ,NZONES,CK    ,CL    ,CKP   ,R1    ,   PART0016
     XZ1    ,UG    ,VG    ,TT0   ,DG    ,UGR   ,UGZ   ,VRG   ,VGZ   ,   PART0017
     XZW    ,R234  ,WPWGT ,WPWT  ,CKK   ,CLL   ,USG   ,ZW1   ,THJ   ,   PART0018
     XSAUR                                                              PART0019
        COMMON/ERR/IERR                                                 PART0020
      COMMON/NAMEA/VAR(6),IORDER,SOWCON,ISW1,ISW2                       PART0021
      COMMON/NAME1/THF,MARK,UUS,UDUDS,ZSINK                             PART0022
      COMMON/NAME3/A20(99)                                              PART0023
C                                                                       PART0024
C         GENERATE ZONE COEFICIENTS                                     PART0025
C                                                                       PART0026
      USG=SQRT ((VK-1.)/(VK+1.)*2.*CPG*TG0/BC)                          PART0027
      WRITE (6,900)VK,WDG                                               PART0028
      C1THI=1.-COS (THI)                                                PART0029
      PSICC=.15915494/RT**2/C1THI                                       PART0030
      STHI=SIN (THI)                                                    PART0031
      RBAR=(1.+RRT*C1THI)/STHI                                          PART0032
      ZSINK=-RRT*STHI-RBAR*C1THI                                        PART0033
      IZ=ZI                                                             PART0034
      JZ=ZJ                                                             PART0035
      NZONES=IZ+JZ+1                                                    PART0036
C                                                                       PART0037
      DO 23 JJ=1,3                                                      PART0038
C                                                                       PART0039
      CALL ONED                                                         PART0040
      IF(IZ) 2,1,2                                                      PART0041
 1    TH=0.                                                             PART0042
      GO TO 201                                                         PART0043
 2    DTH1=THI/ZI                                                       PART0044
      TH=-THI                                                           PART0045
      IF(THI-THF) 201,2001,2001                                         PART0046
 2001 DO 2002 J=1,6                                                     PART0047
      SAV(J,1)=VAR(J)                                                   PART0048
 2002 VAR(J)=SAUR(J)                                                    PART0049
      SOWCON=-1.                                                        PART0050
      CALL JAMES(1.,0.,T,VK)                                            PART0051
      IF(IERR .NE. 0) RETURN                                            PART0052
      SOWCON=0.                                                         PART0053
      DO 2003 J=1,6                                                     PART0054
 2003 VAR(J)=SAV(J,1)                                                   PART0055
 201  IF(JZ) 3,301,3                                                    PART0056
 3    DTH2=THJ/ZJ                                                       PART0057
 301  DO 9 I=1,NZONES                                                   PART0058
      RW=1.+RRT*(1.-COS (TH))                                           PART0059
      ZW1(I)=RRT*SIN (TH)                                               PART0060
      CALL JAMES(RW,ZW1(I),T,VK)                                        PART0061
      IF(IERR .NE. 0) RETURN                                            PART0062
      IF(VAR(2)-.1) 305,305,306                                         PART0063
 305  WRITE (6,94)                                                      PART0064
 94   FORMAT(19HERROR....PARTIL,305)                                    PART0065
      IERR=1                                                            PART0066
      RETURN                                                            PART0067
 306  DO 4 J=1,6                                                        PART0068
 4    SAV(J,I)=VAR(J)                                                   PART0069
      DO 5 J=1,53                                                       PART0070
 5    SAVE(J,I)=A20(J)                                                  PART0071
      IF(I-IZ)6,7,8                                                     PART0072
 6    TH=TH+DTH1                                                        PART0073
      GO TO 9                                                           PART0074
 7    TH=0.                                                             PART0075
      GO TO 9                                                           PART0076
 8    TH=TH+DTH2                                                        PART0077
 9    CONTINUE                                                          PART0078
C                                                                       PART0079
C         INTEGRATION ALONG AXIS                                        PART0080
C                                                                       PART0081
      ASSIGN 144 TO NNN                                                 PART0082
      IF(NZONES-1) 904,904,901                                          PART0083
 901  ASSIGN 140 TO NNN                                                 PART0084
      NN=1                                                              PART0085
      DO 902 J=1,6                                                      PART0086
 902  VAR(J)=SAV(J,1)                                                   PART0087
      ZW=ZW1(1)                                                         PART0088
 904  LAST=0                                                            PART0089
      DZ1=DZI                                                           PART0090
      R1=0.                                                             PART0091
      Z1=ZSINK                                                          PART0092
      CALL PROP(1)                                                      PART0093
      DO 10 I=1,ITOT                                                    PART0094
      TP(I)=TG0*(1.-CL(I)*(1.-TT0))                                     PART0095
 10   UP(I)=CK(I)*UG                                                    PART0096
      DO 15 II=1,10000                                                  PART0097
      IF(ZAX-Z1-DZ1) 11,11,12                                           PART0098
 11   LAST=1                                                            PART0099
      DZ1=ZAX-Z1                                                        PART0100
 12   CALL PROP(1)                                                      PART0101
      TG=TG0*TT0                                                        PART0102
      TT0N=TT0**CAPN                                                    PART0103
      Q1=GMG0*TT0N/DG                                                   PART0104
      Q2=SQRT (RCAP*TG)                                                 PART0105
      DO 14 I=1,ITOT                                                    PART0106
      CBAR=Q1/R(I)/Q2                                                   PART0107
      RE=2.*(UG-UP(I))/Q1*R(I)                                          PART0108
      CALL FIND(IFG,FA,RE,FP)                                           PART0109
      F=FP*(1.+2.52*CBAR)/(1.+3.78*CBAR)                                PART0110
      G=FP/(1.+CG*CBAR*FP)                                              PART0111
      A0=A0C*TT0N/R(I)**2*F                                             PART0112
      B0=A0/1.5*CTOP*G/F                                                PART0113
      DD=DZ1/UP(I)                                                      PART0114
      DUP=A0*(UG-UP(I))*DD                                              PART0115
      DTP=-B0/CPL*(TP(I)-TG)*DD                                         PART0116
      UP(I)=UP(I)+DUP                                                   PART0117
 14   TP(I)=TP(I)+DTP                                                   PART0118
      Z1=Z1+DZ1                                                         PART0119
      GO TO NNN, (140,144)                                              PART0120
 140  IF(UG/USG-(VAR(2)+SAV(2,NN+1))/2.) 144,144,141                    PART0121
 141  NN=NN+1                                                           PART0122
      DO 142 J=1,6                                                      PART0123
 142  VAR(J)=SAV(J,NN)                                                  PART0124
      ZW=ZW1(NN)                                                        PART0125
      IF(NN+1-NZONES) 144,144,143                                       PART0126
 143  ASSIGN 144 TO NNN                                                 PART0127
 144  IF(LAST) 16,15,16                                                 PART0128
 15   CONTINUE                                                          PART0129
 16   TG0TG=TG0-TG                                                      PART0130
      DO 17 I=1,ITOT                                                    PART0131
      CKU0(I)=UP(I)/UG                                                  PART0132
 17   CLT0(I)=(TG0-TP(I))/TG0TG                                         PART0133
C                                                                       PART0134
C         CALCULATION OF K BAR=VK                                       PART0135
C                                                                       PART0136
 20   WL=0.                                                             PART0137
      WK=0.                                                             PART0138
      DO 21 I=1,ITOT                                                    PART0139
      WL=WL+WPWG(I)*CLT0(I)                                             PART0140
 21   WK=WK+WPWG(I)*CKU0(I)**2                                          PART0141
      C1=(1.+CPL/CPG*WL)/(1.+WK)                                        PART0142
      BC=1./C1                                                          PART0143
      C2=0.                                                             PART0144
      DO 22 I=1,ITOT                                                    PART0145
 22   C2=C2+WPWG(I)*(CPL/CPG*CLT0(I)+CKU0(I)*(1.-CKU0(I))*C1)           PART0146
      VK=1.+(GAMMA-1.)/(1.+GAMMA*C2)                                    PART0147
      USG=SQRT ((VK-1.)/(VK+1.)*2.*CPG*TG0/BC)                          PART0148
      WRITE (6,90)VK                                                    PART0149
      ZW=ZW1(NZONES)                                                    PART0150
      CALL JAMES(RW,ZW,T,VK)                                            PART0151
      IF(IERR .NE. 0) RETURN                                            PART0152
      CALL WDGI                                                         PART0153
      DO 225 J=1,6                                                      PART0154
 225  VAR(J)=SAV(J,1)                                                   PART0155
 23   CONTINUE                                                          PART0156
C                                                                       PART0157
C         THE PRINCIPLE INTEGRATION LOOP                                PART0158
C                                                                       PART0159
 24   CN=4.                                                             PART0160
      THN=THI/20.                                                       PART0161
      DO 27 J=1,4                                                       PART0162
      COSTH1=1.-COS (THN)                                               PART0163
      R1=RBAR*SIN (THN)                                                 PART0164
      Z1=ZSINK+RBAR*COSTH1                                              PART0165
      CALL TRACE                                                        PART0166
      PSIC=PSICC*COSTH1                                                 PART0167
      DO 26 I=1,ITOT                                                    PART0168
      PSI(I,J)=WPWG(I)*WDG*PSIC                                         PART0169
      CKU(I,J)=CKK(I)                                                   PART0170
      CKV(I,J)=CKP(I)                                                   PART0171
      CLT(I,J)=CLL(I)                                                   PART0172
 26   R2(I,J)=R234(I)**2                                                PART0173
      CN=CN-1.                                                          PART0174
 27   THN=THI/CN                                                        PART0175
C                                                                       PART0176
C         SELECT INITIAL LINE AND INTERPOLATE FOR VALUES                PART0177
C                                                                       PART0178
      DO 28 I=1,ITOT                                                    PART0179
      T1(I)=R2(I,4)-R2(I,3)                                             PART0180
      T2(I)=R2(I,4)-R2(I,2)                                             PART0181
      T3(I)=R2(I,3)-R2(I,2)                                             PART0182
      T4(I)=R2(I,2)+R2(I,3)                                             PART0183
      T5(I)=R2(I,2)*T3(I)*T2(I)                                         PART0184
      T6(I)=R2(I,3)*T3(I)                                               PART0185
 28   T7(I)=R2(I,4)*T2(I)                                               PART0186
      C=(ZAX-PW(2))/PW**2                                               PART0187
      J=0                                                               PART0188
      CM=NILP                                                           PART0189
      CN=CM                                                             PART0190
      RD=PW/CM/2.                                                       PART0191
      K=0                                                               PART0192
      SKIP=0.                                                           PART0193
      ASSIGN 54 TO JWALL                                                PART0194
      NILP1=NILP+1                                                      PART0195
      DO 70 KK=1,NILP1                                                  PART0196
      RL=PW*CN/CM                                                       PART0197
 53   K=K+1                                                             PART0198
      GO TO JWALL,(54,55)                                               PART0199
 54   ASSIGN 55 TO JWALL                                                PART0200
      GO TO 57                                                          PART0201
 55   IF(J+1-ITOT) 551,551,56                                           PART0202
 551  IF(R234(J+1)-RL+RD) 56,56,60                                      PART0203
 56   IF(SKIP) 69,57,69                                                 PART0204
 57   PS(5,K)=RL                                                        PART0205
      PS(6,K)=ZAX-C*RL**2                                               PART0206
      GO TO 61                                                          PART0207
 60   J=J+1                                                             PART0208
      PS(5,K)=R234(J)                                                   PART0209
      PS(6,K)=ZAX-C*R2(J,4)                                             PART0210
      PS(8,K)=J                                                         PART0211
      SKIP=1.                                                           PART0212
 61   R1=PS(5,K)                                                        PART0213
      Z1=PS(6,K)                                                        PART0214
      CALL PROP(1)                                                      PART0215
      PS(1,K)=DG*RCAP*TT0*TG0                                           PART0216
      PS(2,K)=DG                                                        PART0217
      PS(3,K)=UG                                                        PART0218
      PS(4,K)=VG                                                        PART0219
      IF(J) 70,70,62                                                    PART0220
 62   RR=PS(5,K)**2                                                     PART0221
      DO 63 I=1,J                                                       PART0222
      TC=(T1(I)*(CKU(I,2)-CKU0(I))/T5(I)-(CKU(I,3)-CKU0(I))/T6(I)       PART0223
     X   +(CKU(I,4)-CKU0(I))/T7(I))/T1(I)                               PART0224
      TB=((CKU(I,3)-CKU0(I))/R2(I,3)-(CKU(I,2)-CKU0(I))/R2(I,2))        PART0225
     X   /T3(I)-TC*T4(I)                                                PART0226
      TA=(CKU(I,2)-CKU0(I))/R2(I,2)-(TB+TC*R2(I,2))*R2(I,2)             PART0227
      PS(I+30,K)=UG*(CKU0(I)+RR*(TA+RR*(TB+RR*TC)))                     PART0228
      TC=(T1(I)*(CKV(I,2)-CKV(I,1))/T5(I)-(CKV(I,3)-CKV(I,1))/T6(I)     PART0229
     X   +(CKV(I,4)-CKV(I,1))/T7(I))/T1(I)                              PART0230
      TB=((CKV(I,3)-CKV(I,1))/R2(I,3)-(CKV(I,2)-CKV(I,1))/R2(I,2))      PART0231
     X   /T3(I)-TC*T4(I)                                                PART0232
      TA=(CKV(I,2)-CKV(I,1))/R2(I,2)-(TB+TC*R2(I,2))*R2(I,2)            PART0233
      PS(I+40,K)=VG*(CKV(I,1)+RR*(TA+RR*(TB+RR*TC)))                    PART0234
      TC=(T1(I)*(CLT(I,2)-CLT0(I))/T5(I)-(CLT(I,3)-CLT0(I))/T6(I)       PART0235
     X   +(CLT(I,4)-CLT0(I))/T7(I))/T1(I)                               PART0236
      TB=((CLT(I,3)-CLT0(I))/R2(I,3)-(CLT(I,2)-CLT0(I))/R2(I,2))        PART0237
     X   /T3(I)-TC*T4(I)                                                PART0238
      TA=(CLT(I,2)-CLT0(I))/R2(I,2)-(TB+TC*R2(I,2))*R2(I,2)             PART0239
      TPI=TG0+TG0*(TT0-1.)*(CLT0(I)+RR*(TA+RR*(TB+RR*TC)))              PART0240
      PS(I+10,K)=HPL+CPL*(TPI-TPM)                                      PART0241
      TC=(T1(I)*(PSI(I,2)/R2(I,2)-PSI(I,1)/R2(I,1))/T5(I)-(PSI(I,3)/R2(IPART0242
     X   ,3)-PSI(I,1)/R2(I,1))/T6(I)+(PSI(I,4)/R2(I,4)-PSI(I,1)/R2(I,1))PART0243
     X   /T7(I))/T1(I)                                                  PART0244
      TB=((PSI(I,3)/R2(I,3)-PSI(I,1)/R2(I,1))/R2(I,3)-(PSI(I,2)/R2(I,2) PART0245
     X   -PSI(I,1)/R2(I,1))/R2(I,2))/T3(I)-TC*T4(I)                     PART0246
      TA=(PSI(I,2)/R2(I,2)-PSI(I,1)/R2(I,1))/R2(I,2)-(TB+TC*R2(I,2))    PART0247
     X   *R2(I,2)                                                       PART0248
      DRDZ=-2.*C*PS(5,K)                                                PART0249
 63   PS(I+20,K)=(2.*PSI(I,1)/R2(I,1)+RR*(4.*TA+RR*(6.*TB+RR*8.*TC)))   PART0250
     X   /(PS(I+30,K)-PS(I+40,K)*DRDZ)                                  PART0251
      IF(SKIP) 53,70,53                                                 PART0252
 69   SKIP=0.                                                           PART0253
      K=K-1                                                             PART0254
 70   CN=CN-1.                                                          PART0255
      PS(5,K)=0.                                                        PART0256
C                                                                       PART0257
C         PRINT INITIAL LINE                                            PART0258
C                                                                       PART0259
      NIL=K                                                             PART0260
      WRITE (6,91)                                                      PART0261
      J=0                                                               PART0262
      DO 73 I=1,K                                                       PART0263
      WRITE (6,92) (PS(L,I),L=1,6)                                      PART0264
      IF(PS(8,I)) 71,701,71                                             PART0265
 701  IF(J) 73,73,722                                                   PART0266
 71   J=J+1                                                             PART0267
 722  DO 72 JJ=1,J                                                      PART0268
      LL=10+JJ                                                          PART0269
      LLL=40+JJ                                                         PART0270
 72   WRITE (6,93)JJ,(PS(L,I),L=LL,LLL,10)                              PART0271
 73   CONTINUE                                                          PART0272
      DO 74 J=1,6                                                       PART0273
 74   VAR(J)=SAV(J,1)                                                   PART0274
      RETURN                                                            PART0275
C                                                                       PART0276
 900  FORMAT (7H K BAR=F9.7,6X4HWDG=F10.6,10X                           PART0277
     X24H EQUILIBRIUM CONDITIONS. )                                     PART0278
 90   FORMAT (7H K BAR=F9.7)                                            PART0279
 91   FORMAT(1H1,7X,1HP,16X,3HRHO,14X,1HU,16X,1HV,16X,1HR,16X,1HZ,/     PART0280
     X 15X,4HH-PJ,11X,6HRHO-PJ,9X,4HU-PJ,11X,4HV-PJ)                    PART0281
 92   FORMAT (/6E17.7)                                                  PART0282
 93   FORMAT (5X2HJ=I2,4E15.6)                                          PART0283
       END                                                              PART0284
$IBFTC PCALC                                                            PCAL0000
CPCALC                                                                  PCAL0001
      SUBROUTINE PCALC                                                  PCAL0002
C                                                                       PCAL0003
      COMMON/NAMEA/VAR(6),IORDER,SOWCON,ISW1,ISW2                       PCAL0004
      COMMON/NAME1/THF,MARK,UUS,UDUDS,ZSINK                             PCAL0005
      COMMON/NAME2/VK,NFV,NVAR(6),X,U                                   PCAL0006
     1  ,ALPHA,BETA,GAMMA,DELTA,F(6),X2,X3,W(6)                         PCAL0007
     2  ,MF(6),D(6),AB(6,7),BLOCK(200),MAXIT,NUMIT                      PCAL0008
     3  ,X4,NSAV,XW,XW2,XW3,VKK,T2,XSTAR,Y                              PCAL0009
      DIMENSION A20(1),B21(1),C22(1),D23(1)                             PCAL0010
      COMMON/NAME3/                                                     PCAL0011
     1  A20,A21,A40,B21,B22,B40,B41,B42,B60,B61,B80,C22,C23,C40,C41,C42,PCAL0012
     2C43,C60,C61,C62,C63,C80,C81,C82,C100,C101,C120,D23,D24,D41,D42,D43PCAL0013
     3,D44,D60,D61,D62,D63,D64,D80,D81,D82,D83,D84,D100,D101,D102,D103  PCAL0014
     4 ,D120,D121,D122,D140,D141,D160,D40                               PCAL0015
     5 ,R,R2,R3,R4,R5,R6,R7,R8,R9,R10,R11,R12,R13,R14,R15               PCAL0016
     6 ,D1,D2,D3,D4,D5,D6,D7,D8,D9,D10,D11,D12,D13,D14,D15              PCAL0017
     7 ,UF,VF,DRDX1,DRDX2,DRDX3,DRDX4,DRDX5,DUDX1                       PCAL0018
     8 ,DUDX2,DUDX3,DUDX4,DVDX1,DVDX2,DVDX3,DVDX4                       PCAL0019
C                                                                       PCAL0020
      FTX=1.0/(T2-XW2)                                                  PCAL0021
      FTX1=SQRT(FTX)                                                    PCAL0022
      FTX3=FTX*FTX1                                                     PCAL0023
      FTX5=FTX*FTX3                                                     PCAL0024
      FTX7=FTX*FTX5                                                     PCAL0025
      FTX9=FTX*FTX7                                                     PCAL0026
C                                                                       PCAL0027
      E2=2.0*X                                                          PCAL0028
      E3=3.0*X2                                                         PCAL0029
      E4=4.0*X3                                                         PCAL0030
C                                                                       PCAL0031
      AVR2=A20+A21*X                                                    PCAL0032
C                                                                       PCAL0033
      BVR2=B21*X+B22*X2                                                 PCAL0034
      BVR4=B40+B41*X+B42*X2                                             PCAL0035
      BVR6=B60+B61*X                                                    PCAL0036
C                                                                       PCAL0037
      CVR2=C22*X2+C23*X3                                                PCAL0038
      CVR4=C40+C41*X+C42*X2+C43*X3                                      PCAL0039
      CVR6=C60+C61*X+C62*X2+C63*X3                                      PCAL0040
      CVR8=C80+C81*X+C82*X2                                             PCAL0041
      CVR10=C100+C101*X                                                 PCAL0042
C                                                                       PCAL0043
      DVR2=D23*X3+D24*X4                                                PCAL0044
      DVR4=D41*X+D42*X2+D43*X3+D44*X4                                   PCAL0045
      DVR6=D60+D61*X+D62*X2+D63*X3+D64*X4                               PCAL0046
      DVR8=D80+D81*X+D82*X2+D83*X3+D84*X4                               PCAL0047
      DVR10=D100+D101*X+D102*X2+D103*X3                                 PCAL0048
      DVR12=D120+D121*X+D122*X2                                         PCAL0049
      DVR14=D140+D141*X                                                 PCAL0050
C                                                                       PCAL0051
      BUR2=B21+B22*E2                                                   PCAL0052
      BUR4=B41+B42*E2                                                   PCAL0053
C                                                                       PCAL0054
      CUR2=C22*E2+C23*E3                                                PCAL0055
      CUR4=C41+C42*E2+C43*E3                                            PCAL0056
      CUR6=C61+C62*E2+C63*E3                                            PCAL0057
      CUR8=C81+C82*E2                                                   PCAL0058
C                                                                       PCAL0059
      DUR2=D23*E3+D24*E4                                                PCAL0060
      DUR4=D41+D42*E2+D43*E3+D44*E4                                     PCAL0061
      DUR6=D61+D62*E2+D63*E3+D64*E4                                     PCAL0062
      DUR8=D81+D82*E2+D83*E3+D84*E4                                     PCAL0063
      DUR10=D101+D102*E2+D103*E3                                        PCAL0064
      DUR12=D121+D122*E2                                                PCAL0065
C                                                                       PCAL0066
      CUX2=2.0*C22+3.0*E2*C23                                           PCAL0067
      CUX4=2.0*C42+3.0*E2*C43                                           PCAL0068
      CUX6=2.0*C62+3.0*E2*C63                                           PCAL0069
C                                                                       PCAL0070
      DUX2=3.0*E2*D23+4.0*E3*D24                                        PCAL0071
      DUX4=2.0*D42+3.0*E2*D43+4.0*E3*D44                                PCAL0072
      DUX6=2.0*D62+3.0*E2*D63+4.0*E3*D64                                PCAL0073
      DUX8=2.0*D82+3.0*E2*D83+4.0*E3*D84                                PCAL0074
      DUX10=2.0*D102+3.0*E2*D103                                        PCAL0075
C                                                                       PCAL0076
      DXX2=6.0*(D23+2.0*E2*D24)                                         PCAL0077
      DXX4=6.0*(D43+2.0*E2*D44)                                         PCAL0078
      DXX6=6.0*(D63+2.0*E2*D64)                                         PCAL0079
      DXX8=6.0*(D83+2.0*E2*D84)                                         PCAL0080
C                                                                       PCAL0081
      F2=A21+BUR2+CUR2+DUR2                                             PCAL0082
      F4=BUR4+CUR4+DUR4                                                 PCAL0083
      F6=B61+CUR6+DUR6                                                  PCAL0084
      F8=CUR8+DUR8                                                      PCAL0085
      F10=C101+DUR10                                                    PCAL0086
C                                                                       PCAL0087
      G2=AVR2+BVR2+CVR2+DVR2                                            PCAL0088
      G4=A40+BVR4+CVR4+DVR4                                             PCAL0089
      G6=BVR6+CVR6+DVR6                                                 PCAL0090
      G8=B80+CVR8+DVR8                                                  PCAL0091
      G10=CVR10+DVR10                                                   PCAL0092
      G12=C120+DVR12                                                    PCAL0093
C                                                                       PCAL0094
      H2=2.0*B22+CUX2+DUX2                                              PCAL0095
      H4=2.0*B42+CUX4+DUX4                                              PCAL0096
      H6=CUX6+DUX6                                                      PCAL0097
      H8=2.0*C82+DUX8                                                   PCAL0098
C                                                                       PCAL0099
      P2=6.0*C23+DXX2                                                   PCAL0100
      P4=6.0*C43+DXX4                                                   PCAL0101
      P6=6.0*C63+DXX6                                                   PCAL0102
C                                                                       PCAL0103
      UF1R=D1*F2+D3*F4+D5*F6+D7*F8+D9*F10+D11*DUR12+                    PCAL0104
     1D13*D141                                                          PCAL0105
      VF1R=2.0*(G2+2.0*D2*G4+3.0*D4*G6+4.0*D6*G8+5.0*D8*G10+            PCAL0106
     1     6.0*D10*G12+7.0*D12*DVR14+8.0*D14*D160)                      PCAL0107
      UF2R=2.0*(F2+2.0*D2*F4+3.0*D4*F6+4.0*D6*F8+5.0*D8*F10+            PCAL0108
     1     6.0*D10*DUR12+7.0*D12*D141)                                  PCAL0109
      VF2R=2.0*(6.0*D1*G4+15.0*D3*G6+28.0*D5*G8+45.0*D7*G10+            PCAL0110
     1     66.0*D9*G12+91.0*D11*DVR14+120.0*D13*D160)                   PCAL0111
      UF3R=2.0*(6.0*D1*F4+15.0*D3*F6+28.0*D5*F8+45.0*D7*F10+            PCAL0112
     1     66.0*D9*DUR12+91.0*D11*D141)                                 PCAL0113
      VF3R=4.0*(6.0*G4+30.0*D2*G6+84.0*D4*G8+180.0*D6*G10+              PCAL0114
     1     330.0*D8*G12+546.0*D10*DVR14+840.0*D12*D160)                 PCAL0115
      UF4R=4.0*(6.0*F4+30.0*D2*F6+84.0*D4*F8+180.0*D6*F10+              PCAL0116
     1     330.0*D8*DUR12+546.0*D10*D141)                               PCAL0117
      VF4R=8.0*(45.0*D1*G6+210.0*D3*G8+630.0*D5*G10+1485.0*D7*G12+      PCAL0118
     1      3003.0*D9*DVR14+5460.0*D11*D160)                            PCAL0119
      UF1X=ALPHA+BETA*X+(GAMMA*X2)/2.0+(DELTA*X3)/6.0+R2*H2+R4*H4+      PCAL0120
     1   R6*H6+R8*H8+R10*DUX10+2.0*R12*D122                             PCAL0121
      VF1X=D1*F2+D3*F4+D5*F6+D7*F8+D9*F10+D11*DUR12+D13*D141            PCAL0122
      UF1X1R=D1*H2+D3*H4+D5*H6+D7*H8+D9*DUX10+2.0*D11*D122              PCAL0123
      VF1X1R=2.0*(F2+2.0*D2*F4+3.0*D4*F6+4.0*D6*F8+5.0*D8*F10+6.0*D10   PCAL0124
     1      *DUR12+7.0*D12*D141)                                        PCAL0125
      UF1X2R=2.0*(H2+2.0*D2*H4+3.0*D4*H6+4.0*D6*H8+5.0*D8*DUX10+12.0*   PCAL0126
     1      D10*D122)                                                   PCAL0127
      VF1X2R=2.0*(6.0*D1*F4+15.0*D3*F6+28.0*D5*F8+45.0*D7*F10+66.0*D9*  PCAL0128
     1      DUR12+91.0*D11*D141)                                        PCAL0129
      UF1X3R=2.0*(6.0*D1*H4+15.0*D3*H6+28.0*D5*H8+45.0*D7*DUX10+132.0*  PCAL0130
     1      D9*D122)                                                    PCAL0131
      VF1X3R=12.0*(2.0*F4+10.0*D2*F6+28.0*D4*F8+60.0*D6*F10+110.0*      PCAL0132
     1      D8*DUR12+182.0*D10*D141)                                    PCAL0133
      UF2X=BETA+GAMMA*X+(DELTA*X2)/2.0+R2*P2+R4*P4+R6*P6+               PCAL0134
     1 R8*DXX8+6.0*R10*D103                                             PCAL0135
      VF2X=D1*H2+D3*H4+D5*H6+D7*H8+D9*DUX10+2.0*D11*D122                PCAL0136
      UF2X1R=D1*P2+D3*P4+D5*P6+D7*DXX8+6.0*D9*D103                      PCAL0137
      VF2X1R=2.0*(H2+2.0*D2*H4+3.0*D4*H6+4.0*D6*H8+5.0*D8*DUX10+        PCAL0138
     1     12.0*D10*D122)                                               PCAL0139
      UF2X2R=2.0*(P2+2.0*D2*P4+3.0*D4*P6+4.0*D6*DXX8+30.0*D8*D103)      PCAL0140
      VF2X2R=2.0*(6.0*D1*H4+15.0*D3*H6+28.0*D5*H8+45.0*D7*DUX10+        PCAL0141
     1     132.0*D9*D122)                                               PCAL0142
      UF3X=GAMMA+DELTA*X+24.0*(R2*D24+R4*D44+R6*D64+R8*D84)             PCAL0143
      VF3X=D1*P2+D3*P4+D5*P6+D7*DXX8+6.0*D9*D103                        PCAL0144
      UF3X1R=24.0*(D1*D24+D3*D44+D5*D64+D7*D84)                         PCAL0145
      VF3X1R=2.0*(P2+2.0*D2*P4+3.0*D4*P6+4.0*D6*DXX8+30.0*              PCAL0146
     1      D8*D103)                                                    PCAL0147
      UF4X=DELTA                                                        PCAL0148
      VF4X=24.0*(D1*D24+D3*D44+D5*D64+D7*D84)                           PCAL0149
C                                                                       PCAL0150
      UF=ALPHA*X+.5*BETA*X2+(GAMMA*X3)/6.0+(DELTA*X4)/24.0+             PCAL0151
     1  R2*F2+R4*F4+R6*F6+R8*F8+R10*F10+R12*DUR12+R14*D141              PCAL0152
      VF=D1*G2+D3*G4+D5*G6+D7*G8+D9*G10+D11*G12+                        PCAL0153
     1D13*DVR14+D15*D160                                                PCAL0154
C                                                                       PCAL0155
      DRDX1=XW*FTX1                                                     PCAL0156
      DRDX2=FTX1+XW2*FTX3                                               PCAL0157
      DRDX3=3.0*XW*FTX3+3.0*XW3*FTX5                                    PCAL0158
      DRDX4=3.0*FTX3+18.0*XW2*FTX5+15.0*XW2*XW2*FTX7                    PCAL0159
      DRDX5=45.0*XW*FTX5+150.0*XW3*FTX7+105.0*XW2*XW3*FTX9              PCAL0160
      DRDXS=DRDX1*DRDX1                                                 PCAL0161
      DUDX1=UF1X+UF1R*DRDX1                                             PCAL0162
      DUDX2=UF2X+2.0*UF1X1R*DRDX1+UF2R*DRDXS+UF1R*DRDX2                 PCAL0163
      DUDX3=UF3X+3.0*UF2X1R*DRDX1+DRDXS*(3.0*UF1X2R+UF3R*DRDX1)+        PCAL0164
     1  3.0*DRDX2*(UF1X1R+UF2R*DRDX1)+UF1R*DRDX3                        PCAL0165
      DUDX4=UF4X+4.0*UF3X1R*DRDX1+3.0*UF2X1R*DRDX2+DRDXS*               PCAL0166
     1  (6.0*UF2X2R+4.0*UF1X3R*DRDX1+UF4R*DRDXS)+DRDX1*DRDX2*           PCAL0167
     2     (12.0*UF1X2R+6.0*UF3R*DRDX1)+DRDX2*(3.0*UF2X1R+3.0*UF2R*DRDX2PCAL0168
     3)+      4.0*DRDX3*(UF1X1R+UF2R*DRDX1)+UF1R*DRDX4                  PCAL0169
      DVDX1=VF1X+VF1R*DRDX1                                             PCAL0170
      DVDX2=VF2X+2.0*VF1X1R*DRDX1+VF2R*DRDXS+VF1R*DRDX2                 PCAL0171
      DVDX3=VF3X+3.0*VF2X1R*DRDX1+DRDXS*(3.0*VF1X2R+VF3R*DRDX1)+        PCAL0172
     1  3.0*DRDX2*(VF1X1R+VF2R*DRDX1)+VF1R*DRDX3                        PCAL0173
      DVDX4=VF4X+4.0*VF3X1R*DRDX1+3.0*VF2X1R*DRDX2+DRDXS*               PCAL0174
     1   (6.0*VF2X2R+4.0*VF1X3R*DRDX1+VF4R*DRDXS)+DRDX1*DRDX2*          PCAL0175
     2     (12.0*VF1X2R+6.0*VF3R*DRDX1)+DRDX2*3.0*(VF2X1R+VF2R*DRDX2)+  PCAL0176
     3       4.0*DRDX3*(VF1X1R+VF2R*DRDX1)+VF1R*DRDX4                   PCAL0177
 999  RETURN                                                            PCAL0178
C                                                                       PCAL0179
      END                                                               PCAL0180
$IBFTC PROP                                                             PROP0000
      SUBROUTINE PROP(I)                                                PROP0001
C                                                                       PROP0002
      DIMENSION  WPWG(10),PS(50,50),PW(150),R(10),XXXX(4881),FA(100),   PROP0003
     XGF(10),WDP(10),CK(10),CL(10),CKP(10),WPWT(10),SAV(6,10),          PROP0004
     XSAVE(53,10),R234(10),CKK(10),CLL(10)                              PROP0005
      COMMON CAPN  ,CPL   ,CPS   ,CTOP  ,DL    ,DR    ,EK    ,EKW   ,   PROP0006
     XEW    ,WPWG  ,GAMMA ,GMG0  ,HPL   ,HPS   ,IMAX  ,N1    ,N2    ,   PROP0007
     XPS    ,PG0   ,PW    ,R     ,RCAP  ,RHOP0 ,SMP   ,TG0   ,TPM   ,   PROP0008
     XSAV   ,SAVE  ,XXXX  ,FA    ,VGR   ,CPG   ,NILP  ,THI   ,GF    ,   PROP0009
     XWDG   ,IFG   ,RHOG0 ,THRUST,TI    ,WDP   ,TMD   ,IPPX  ,BC    ,   PROP0010
     XTHII  ,DZMIN ,ITOT  ,ZAX   ,RRT   ,RT    ,VK    ,T     ,CG    ,   PROP0011
     XA0C   ,DZI   ,ZI    ,ZJ    ,NZONES,CK    ,CL    ,CKP   ,R1    ,   PROP0012
     XZ1    ,UG    ,VG    ,TT0   ,DG    ,UGR   ,UGZ   ,VRG   ,VGZ   ,   PROP0013
     XZW    ,R234  ,WPWGT ,WPWT  ,CKK   ,CLL   ,USG                     PROP0014
      COMMON/NAMEA/VAR(6),IORDER,SOWCON,ISW1,ISW2                       PROP0015
      COMMON/NAME3/                                                     PROP0016
     X  A20,A21,A40,B21,B22,B40,B41,B42,B60,B61,B80,C22,C23,C40,C41,C42,PROP0017
     X   C43,C60,C61,C62,C63,C80,C81,C82,C100,C101,C120                 PROP0018
     1 ,DUMMY(72)                                                       PROP0019
C                                                                       PROP0020
      E=VK-1.                                                           PROP0021
      Q3=E/(VK+1.)                                                      PROP0022
      Q2=Z1+VAR(1)-ZW                                                   PROP0023
      P0=VAR(2)+Q2*(VAR(3)+Q2*(VAR(4)/2.+Q2*VAR(5)/6.))                 PROP0024
      IF(R1) 20,10,20                                                   PROP0025
C                                                                       PROP0026
C         CALCULATE PROPERTIES WITH R1=0                                PROP0027
C                                                                       PROP0028
 10   GO TO (11,12),I                                                   PROP0029
 11   UG=USG*P0                                                         PROP0030
      VG=0.                                                             PROP0031
      TT0=1.-Q3/USG/USG*UG*UG                                           PROP0032
      DG=RHOG0*TT0**(1./E)                                              PROP0033
      GO TO 999                                                         PROP0034
C                                                                       PROP0035
C          CALCULATE DERIVATIVES WITH R1=0                              PROP0036
C                                                                       PROP0037
 12   G0=VAR(3)+Q2*(VAR(4)+Q2*VAR(5)/2.)                                PROP0038
      UGZ=USG*G0                                                        PROP0039
      P1=A20+Q2*(A21+B21+Q2*(B22+C22+Q2*C23))                           PROP0040
      VRG=2.*USG*P1                                                     PROP0041
      UGR=0.                                                            PROP0042
      VGZ=0.                                                            PROP0043
      GO TO 999                                                         PROP0044
C                                                                       PROP0045
 20   P2=A21+B21+Q2*(2.*(B22+C22)+3.*Q2*C23)                            PROP0046
      P4=B41+C41+Q2*(2.*(B42+C42)+Q2*3.*C43)                            PROP0047
      P6=B61+C61+Q2*(2.*C62+Q2*3.*C63)                                  PROP0048
      P8=C81+2.*Q2*C82                                                  PROP0049
      P1=A20+Q2*(A21+B21+Q2*(B22+C22+Q2*C23))                           PROP0050
      P3=A40+B40+C40+Q2*(B41+C41+Q2*(B42+C42+Q2*C43))                   PROP0051
      P5=B60+C60+Q2*(B61+C61+Q2*(C62+Q2*C63))                           PROP0052
      P7=B80+C80+Q2*(C81+Q2*C82)                                        PROP0053
      P9=C100+Q2*C101                                                   PROP0054
      Q4=R1*R1                                                          PROP0055
      GO TO (21,22),I                                                   PROP0056
C                                                                       PROP0057
C         CALCULATE PROPERTIES                                          PROP0058
C                                                                       PROP0059
 21   UG=USG*(P0+Q4*(P2+Q4*(P4+Q4*(P6+Q4*(P8+Q4*C101)))))               PROP0060
      VG=2.*USG*(R1*(P1+Q4*(2.*P3+Q4*(3.*P5+Q4*(4.*P7+Q4*(5.*P9+Q4*6.   PROP0061
     X   *C120))))))                                                    PROP0062
      TT0=1.-Q3/USG/USG*(UG*UG+VG*VG)                                   PROP0063
      DG=RHOG0*TT0**(1./E)                                              PROP0064
      GO TO 999                                                         PROP0065
C                                                                       PROP0066
C          CALCULATE DERIVATIVES                                        PROP0067
C                                                                       PROP0068
 22   UGR=USG*R1*(2.*P2+Q4*(4.*P4+Q4*(6.*P6+Q4*(8.*P8+10.*Q4*C101))))   PROP0069
      VRG=USG*(2.*P1+Q4*(12.*P3+Q4*(30.*P5+Q4*(56.*P7+Q4*(90.*P9+Q4*    PROP0070
     X   132.*C120)))))                                                 PROP0071
      G0=VAR(3)+Q2*(VAR(4)+Q2*VAR(5)/2.)                                PROP0072
      G2=2.*(B22+C22)+Q2*6.*C23                                         PROP0073
      G4=2.*(B42+C42)+Q2*6.*C43                                         PROP0074
      G6=C62+C62+Q2*6.*C63                                              PROP0075
      G8=C82+C82                                                        PROP0076
      UGZ=USG*(G0+Q4*(G2+Q4*(G4+Q4*(G6+Q4*G8))))                        PROP0077
      VGZ=UGR                                                           PROP0078
C                                                                       PROP0079
 999  RETURN                                                            PROP0080
       END                                                              PROP0081
$IBFTC TRACE                                                            TRAC0000
      SUBROUTINE TRACE                                                  TRAC0001
C                                                                       TRAC0002
      DIMENSION  WPWG(10),PS(50,50),PW(150),R(10),XXXX(4881),FA(100),   TRAC0003
     XGF(10),WDP(10),CK(10),CL(10),CKP(10),WPWT(10),SAV(6,10),          TRAC0004
     XSAVE(53,10),R234(10),CKK(10),CLL(10),ZW1(10)                      TRAC0005
     X,RR(10),Z1P(10),Z2P(10),TP(10),UP(10),VP(10)                      TRAC0006
      COMMON CAPN  ,CPL   ,CPS   ,CTOP  ,DL    ,DR    ,EK    ,EKW   ,   TRAC0007
     XEW    ,WPWG  ,GAMMA ,GMG0  ,HPL   ,HPS   ,IMAX  ,N1    ,N2    ,   TRAC0008
     XPS    ,PG0   ,PW    ,R     ,RCAP  ,RHOP0 ,SMP   ,TG0   ,TPM   ,   TRAC0009
     XSAV   ,SAVE  ,XXXX  ,FA    ,VGR   ,CPG   ,NILP  ,THI   ,GF    ,   TRAC0010
     XWDG   ,IFG   ,RHOG0 ,THRUST,TI    ,WDP   ,TMD   ,IPPX  ,BC    ,   TRAC0011
     XTHII  ,DZMIN ,ITOT  ,ZAX   ,RRT   ,RT    ,VK    ,T     ,CG    ,   TRAC0012
     XA0C   ,DZI   ,ZI    ,ZJ    ,NZONES,CK    ,CL    ,CKP   ,R1    ,   TRAC0013
     XZ1    ,UG    ,VG    ,TT0   ,DG    ,UGR   ,UGZ   ,VRG   ,VGZ   ,   TRAC0014
     XZW    ,R234  ,WPWGT ,WPWT  ,CKK   ,CLL   ,USG   ,ZW1              TRAC0015
      COMMON/NAMEA/VAR(6),IORDER,SOWCON,ISW1,ISW2                       TRAC0016
      COMMON/NAME3/A20(99)                                              TRAC0017
C                                                                       TRAC0018
      ASSIGN 23 TO NNN                                                  TRAC0019
      IF(NZONES-1) 9,9,5                                                TRAC0020
 5    ASSIGN 18 TO NNN                                                  TRAC0021
      L=1                                                               TRAC0022
      DO 6 J=1,6                                                        TRAC0023
 6    VAR(J)=SAV(J,1)                                                   TRAC0024
      DO 7 J=1,53                                                       TRAC0025
 7    A20(J)=SAVE(J,1)                                                  TRAC0026
      ZW=ZW1(1)                                                         TRAC0027
 9    ITEST=(ITOT+1)/2                                                  TRAC0028
      CALL PROP(1)                                                      TRAC0029
      DO 10 I=1,ITOT                                                    TRAC0030
      RR(I)=R1                                                          TRAC0031
      Z2P(I)=Z1                                                         TRAC0032
      TP(I)=TG0*(1.-CL(I)*(1.-TT0))                                     TRAC0033
      VP(I)=CK(I)*VG                                                    TRAC0034
 10   UP(I)=CK(I)*UG                                                    TRAC0035
      C=(ZAX-PW(2))/PW**2                                               TRAC0036
C                                                                       TRAC0037
      M=1                                                               TRAC0038
      DO 24 II=1,10000                                                  TRAC0039
      N=M                                                               TRAC0040
      DO 23 I=N,ITOT                                                    TRAC0041
      R1=RR(I)                                                          TRAC0042
      CALL PROP(1)                                                      TRAC0043
      TG=TG0*TT0                                                        TRAC0044
      TT0N=TT0**CAPN                                                    TRAC0045
      Q1=GMG0*TT0N/DG                                                   TRAC0046
      Q2=SQRT (RCAP*TG)                                                 TRAC0047
      CBAR=Q1/R(I)/Q2                                                   TRAC0048
      RE=2.*SQRT ((UG-UP(I))**2+(VG-VP(I))**2)/Q1*R(I)                  TRAC0049
      CALL FIND(IFG,FA,RE,FP)                                           TRAC0050
      F=FP*(1.+2.52*CBAR)/(1.+3.78*CBAR)                                TRAC0051
      G=FP/(1.+CG*CBAR*FP)                                              TRAC0052
      A0=A0C*TT0N/R(I)**2*F                                             TRAC0053
      B0=A0/1.5*CTOP*G/F                                                TRAC0054
      DZ1=DZI                                                           TRAC0055
      DD=DZ1/UP(I)                                                      TRAC0056
      DR1=VP(I)*DD                                                      TRAC0057
      Z1P(I)=Z2P(I)                                                     TRAC0058
      Z2P(I)=ZAX-C*(R1+DR1)**2                                          TRAC0059
      IF(Z1+DZ1-Z2P(I)) 14,14,13                                        TRAC0060
 13   ZP=(Z1-Z1P(I)*(DZ1/(Z2P(I)-Z1P(I))))/(1.-(DZ1/(Z2P(I)-Z1P(I))))   TRAC0061
      DZ1=ZP-Z1                                                         TRAC0062
      DD=DZ1/UP(I)                                                      TRAC0063
      DR1=VP(I)*DD                                                      TRAC0064
      R234(I)=R1+DR1                                                    TRAC0065
      M=I+1                                                             TRAC0066
 14   DUP=A0*(UG-UP(I))*DD                                              TRAC0067
      DVP=A0*(VG-VP(I))*DD                                              TRAC0068
      DTP=-B0/CPL*(TP(I)-TG)*DD                                         TRAC0069
      RR(I)=RR(I)+DR1                                                   TRAC0070
      UP(I)=UP(I)+DUP                                                   TRAC0071
      VP(I)=VP(I)+DVP                                                   TRAC0072
      TP(I)=TP(I)+DTP                                                   TRAC0073
      CKK(I)=UP(I)/UG                                                   TRAC0074
      CKP(I)=VP(I)/VG                                                   TRAC0075
      CLL(I)=(TG0-TP(I))/(TG0-TG)                                       TRAC0076
      IF(I-ITEST) 23,17,23                                              TRAC0077
 17   GO TO NNN,(18,23)                                                 TRAC0078
 18   IF(UG/USG-(VAR(2)+SAV(2,L+1))/2.) 23,23,19                        TRAC0079
 19   L=L+1                                                             TRAC0080
      DO 20 J=1,6                                                       TRAC0081
 20   VAR(J)=SAV(J,L)                                                   TRAC0082
      DO 21 J=1,53                                                      TRAC0083
 21   A20(J)=SAVE(J,L)                                                  TRAC0084
      ZW=ZW1(L)                                                         TRAC0085
      IF(L+1-NZONES) 23,23,22                                           TRAC0086
 22   ASSIGN 23 TO NNN                                                  TRAC0087
 23   CONTINUE                                                          TRAC0088
      Z1=Z1+DZI                                                         TRAC0089
      IF(M-ITOT) 24,24,25                                               TRAC0090
 24   CONTINUE                                                          TRAC0091
 25   RETURN                                                            TRAC0092
       END                                                              TRAC0093
$IBFTC WDGI                                                             WDGI0000
      SUBROUTINE WDGI                                                   WDGI0001
C                                                                       WDGI0002
      DIMENSION  WPWG(10),PS(50,50),PW(150),R(10),XXXX(4881),FA(100),   WDGI0003
     XGF(10),WDP(10),CK(10),CL(10),CKP(10),WPWT(10),SAV(6,10),          WDGI0004
     XSAVE(53,10),R234(10)                                              WDGI0005
      COMMON CAPN  ,CPL   ,CPS   ,CTOP  ,DL    ,DR    ,EK    ,EKW   ,   WDGI0006
     XEW    ,WPWG  ,GAMMA ,GMG0  ,HPL   ,HPS   ,IMAX  ,N1    ,N2    ,   WDGI0007
     XPS    ,PG0   ,PW    ,R     ,RCAP  ,RHOP0 ,SMP   ,TG0   ,TPM   ,   WDGI0008
     XSAV   ,SAVE  ,XXXX  ,FA    ,VGR   ,CPG   ,NILP  ,THI   ,GF    ,   WDGI0009
     XWDG   ,IFG   ,RHOG0 ,THRUST,TI    ,WDP   ,TMD   ,IPPX  ,BC    ,   WDGI0010
     XTHII  ,DZMIN ,ITOT  ,ZAX   ,RRT   ,RT    ,VK    ,T     ,CG    ,   WDGI0011
     XA0C   ,DZI   ,ZI    ,ZJ    ,NZONES,CK    ,CL    ,CKP   ,R1    ,   WDGI0012
     XZ1    ,UG    ,VG    ,TT0   ,DG    ,UGR   ,UGZ   ,VRG   ,VGZ   ,   WDGI0013
     XZW    ,R234  ,WPWGT ,WPWT                                         WDGI0014
C                                                                       WDGI0015
      WDG=0.                                                            WDGI0016
      C1=(ZAX-PW(2))/PW(1)**2                                           WDGI0017
      R1=PW(1)                                                          WDGI0018
      Z1=PW(2)                                                          WDGI0019
      CALL PROP(1)                                                      WDGI0020
      R0=R1                                                             WDGI0021
      Z0=Z1                                                             WDGI0022
      U0=UG                                                             WDGI0023
      V0=VG                                                             WDGI0024
      DG0=DG                                                            WDGI0025
      ASSIGN 3 TO NEND                                                  WDGI0026
      DO 5 I=1,500                                                      WDGI0027
      R1=R0-.02                                                         WDGI0028
      IF(R1) 1,1,2                                                      WDGI0029
 1    R1=0.                                                             WDGI0030
      ASSIGN 6 TO NEND                                                  WDGI0031
 2    Z1=ZAX-C1*R1**2                                                   WDGI0032
      CALL PROP(1)                                                      WDGI0033
      RR=R1-R0                                                          WDGI0034
      WDG=WDG+(DG0*U0+DG*UG-(Z1-Z0)/RR*(DG0*V0+DG*VG))*(R0+R1)*RR       WDGI0035
      GO TO NEND,(3,6)                                                  WDGI0036
 3    R0=R1                                                             WDGI0037
      Z0=Z1                                                             WDGI0038
      U0=UG                                                             WDGI0039
      V0=VG                                                             WDGI0040
 5    DG0=DG                                                            WDGI0041
 6    WDG=-1.5729633*WDG*RT*RT                                          WDGI0042
      DO 20 J=1,ITOT                                                    WDGI0043
 20   WDP(J)=WPWG(J)*WDG                                                WDGI0044
 22   WRITE (6,99)WDG,(WDP(I),I=1,ITOT)                                 WDGI0045
 99   FORMAT(5H0WDG=F9.5,10H   WDP(I)=(10F9.5))                         WDGI0046
      RETURN                                                            WDGI0047
       END                                                              WDGI0048
$ORIGIN        LEVEL1,UT2
$IBFTC N2MAIN                                                           N2MA0000
      SUBROUTINE N2MAIN                                                 N2MA0001
        COMMON/ERR/IERR                                                 N2MA0002
      COMMON/NAMEP/P(50,100)                                            N2MA0003
C                                                                       N2MA0004
      CALL N3MAIN                                                       N2MA0005
      IF(IERR .NE. 0) RETURN                                            N2MA0006
      CALL CNTRL                                                        N2MA0007
C                                                                       N2MA0008
      RETURN                                                            N2MA0009
       END                                                              N2MA0010
$ORIGIN        LEVEL2,UT2
$IBFTC CONSTS                                                           CONS0000
      SUBROUTINE CONSTS                                                 CONS0001
C                                                                       CONS0002
      COMMON/NAMEP/P(50,100)                                            CONS0003
      DIMENSION   A(10)    ,BS(10)   ,FJ(10)             ,PP(50,100),   CONS0004
     X PS(50,50) ,PTB(50,8),PW(150)  ,R(10)    ,TRL1(5)  ,TRL2(5)   ,   CONS0005
     X TRL3(5)   ,TRL4(5)  ,TRM1(5)  ,TRM2(5)  ,TRM3(5)  ,TRM4(5)   ,   CONS0006
     X FA(100)   ,WPWG(10) ,WDP(10)  ,AA(8)                             CONS0007
      COMMON CAPN  ,CPL   ,CPS   ,CTOP  ,DL    ,DR    ,EK    ,EKW   ,   CONS0008
     XEW    ,FJ    ,GAMMA ,GMG0  ,HPL   ,HPS   ,IMAX  ,N1    ,N2    ,   CONS0009
     XPS    ,PG0   ,PW    ,R     ,RCAP  ,RHOP0 ,SMP   ,TG0   ,TPM   ,   CONS0010
     XA     ,BS    ,CC    ,CONST ,EOW1  ,EOW2  ,HIT   ,ICNT  ,IPP   ,   CONS0011
     XJMAX  ,MW           ,PP    ,PTB   ,TRL1  ,TRL2  ,TRL3  ,TRL4  ,   CONS0012
     XTRM1  ,TRM2  ,TRM3  ,TRM4  ,LRC   ,N3    ,FA    ,VGR   ,CPG   ,   CONS0013
     XNILP  ,THI   ,WPWG  ,WDG   ,RC    ,RHOG0 ,THRUST,TI    ,WDP   ,   CONS0014
     XAA    ,RT    ,VK    ,T                                            CONS0015
C                                                                       CONS0016
      DO 10 I=1,50                                                      CONS0017
      DO 10 J=1,100                                                     CONS0018
      P (I,J)=0.                                                        CONS0019
 10   PP(I,J)=0.                                                        CONS0020
      DO 100 I=1,50                                                     CONS0021
 100  PTB(I,3)=0.                                                       CONS0022
      WRITE (6,90)                                                      CONS0023
 90   FORMAT(1H152X17HGAS-PARTICLE FLOW/27H CHARACTERISTIC CALCULATION) CONS0024
C                                                                       CONS0025
C         INTEGRATION OF INITIAL LINE FOR THRUST                        CONS0026
C                                                                       CONS0027
      JJ=JMAX                                                           CONS0028
      TI=0.                                                             CONS0029
      DO 6  I=1,50                                                      CONS0030
      PS(7,I)=(PS(1,I)/PS(2,I))**CAPN                                   CONS0031
      I=I                                                               CONS0032
      IF(PS(1,I+1)) 1,11,1                                              CONS0033
 1    RRRR=(PS(5,I+1)+PS(5,I))*(PS(5,I+1)-PS(5,I))                      CONS0034
      ZZRR=(PS(6,I+1)-PS(6,I))/(PS(5,I+1)-PS(5,I))                      CONS0035
      TI=TI+(PS(1,I)+PS(1,I+1)+PS(2,I)*PS(3,I)**2+                      CONS0036
     X   PS(2,I+1)*PS(3,I+1)**2-ZZRR*                                   CONS0037
     X   (PS(2,I)*PS(3,I)*PS(4,I)+PS(2,I+1)*PS(3,I+1)*PS(4,I+1)))*RRRR  CONS0038
      IF(JJ) 4,4,2                                                      CONS0039
 2    DO 3 J=1,JJ                                                       CONS0040
 3    TI=TI+(PS(J+20,I)*PS(J+30,I)**2+PS(J+20,I+1)*PS(J+30,I+1)**2      CONS0041
     X   -ZZRR*(PS(J+20,I)*PS(J+30,I)*PS(J+40,I)+PS(J+20,I+1)*          CONS0042
     X   PS(J+30,I+1)*PS(J+40,I+1)))*RRRR                               CONS0043
 4    IF(PS(8,I+1)) 5,6 ,5                                              CONS0044
 5    JJ=JJ+1                                                           CONS0045
 6    CONTINUE                                                          CONS0046
 11   THRUST=-.04888927*TI*RT*RT                                        CONS0047
C                                                                       CONS0048
C         CALCULATION OF V SUB G/R AT INITIAL LINE AXIS POINT           CONS0049
C                                                                       CONS0050
      VGR=(PS(4,I-1)/PS(5,I-1))*PS(3,I)/(PS(3,I-1)+                     CONS0051
     X    (PS(6,I)-PS(6,I-1))*PS(4,I-1)/PS(5,I-1))                      CONS0052
      RETURN                                                            CONS0053
      END                                                               CONS0054
$IBFTC N3MAIN                                                           N3MA0000
      SUBROUTINE N3MAIN                                                 N3MA0001
C                                                                       N3MA0002
      CALL CONSTS                                                       N3MA0003
      CALL WALL                                                         N3MA0004
      RETURN                                                            N3MA0005
C                                                                       N3MA0006
      END                                                               N3MA0007
$IBFTC WALL                                                             WALL0000
      SUBROUTINE WALL                                                   WALL0001
C                                                                       WALL0002
      DIMENSION  WPWG(10),PS(50,50),PW(150),R(10),XXXX(4881),FA(100),   WALL0003
     XGF(10),WDP(10),CK(10),CL(10),CKP(10),WPWT(10),SAV(6,10),          WALL0004
     XSAVE(53,10),R234(10),XXX(31),SAUR(6)                              WALL0005
      COMMON CAPN  ,CPL   ,CPS   ,CTOP  ,DL    ,DR    ,EK    ,EKW   ,   WALL0006
     XEW    ,WPWG  ,GAMMA ,GMG0  ,HPL   ,HPS   ,IMAX  ,N1    ,N2    ,   WALL0007
     XPS    ,PG0   ,PW    ,R     ,RCAP  ,RHOP0 ,SMP   ,TG0   ,TPM   ,   WALL0008
     XSAV   ,SAVE  ,XXXX  ,FA    ,VGR   ,CPG   ,NILP  ,THI   ,GF    ,   WALL0009
     XWDG   ,IFG   ,RHOG0 ,THRUST,TI    ,WDP   ,TMD   ,IPPX  ,BC    ,   WALL0010
     XTHII  ,DZMIN ,ITOT  ,ZAX   ,RRT   ,RT    ,VK    ,T     ,CG    ,   WALL0011
     XA0C   ,DZI   ,ZI    ,ZJ    ,NZONES,CK    ,CL    ,CKP   ,R1    ,   WALL0012
     XZ1    ,UG    ,VG    ,TT0   ,DG    ,UGR   ,UGZ   ,VRG   ,VGZ   ,   WALL0013
     XZW    ,R234  ,WPWGT ,WPWT  ,XXX   ,THJ   ,SAUR                    WALL0014
     X,EPS,RWMAX,ZWMAX                                                  WALL0015
        COMMON/ERR/IERR                                                 WALL0016
      COMMON/NAMEW/IWALL,THIR,THJR                                      WALL0017
C                                                                       WALL0018
      IF(IWALL .EQ. 0) RETURN                                           WALL0019
      IF(IWALL .GT. 3 .OR. IWALL .LT. 1) GO TO 999                      WALL0020
      DLTHR=(THJR-THIR)/29.                                             WALL0021
      ANG=THIR                                                          WALL0022
      DO 9 I=1,57,2                                                     WALL0023
      SN=SIN (ANG)                                                      WALL0024
      CS=COS (ANG)                                                      WALL0025
      PW(I)=1.+RRT*(1.-CS)                                              WALL0026
      PW(I+1)=RRT*SN                                                    WALL0027
 9    ANG=ANG+DLTHR                                                     WALL0028
      SN=SIN (THJR)                                                     WALL0029
      CS=COS (THJR)                                                     WALL0030
      PW(59)=1.+RRT*(1.-CS)                                             WALL0031
      PW(60)=RRT*SN                                                     WALL0032
      GO TO(1,2,3),IWALL                                                WALL0033
C                                                                       WALL0034
C         CONE     (IWALL=1)                                            WALL0035
C                                                                       WALL0036
 1    PW(61)=SQRT (EPS)                                                 WALL0037
      PW(62)=(PW(61)-PW(59))*CS/SN+PW(60)                               WALL0038
      J=62                                                              WALL0039
      GO TO 4                                                           WALL0040
C                                                                       WALL0041
C         PARABOLA (IWALL=2)                                            WALL0042
C                                                                       WALL0043
 2    CN1=SN/CS                                                         WALL0044
      CN4=CN1*(ZWMAX-PW(60))                                            WALL0045
      A1=(RWMAX**2-PW(59)**2-2.*PW(59)*CN4)/2./                         WALL0046
     X    (RWMAX-PW(59)-CN4)                                            WALL0047
      B1=2.*(PW(59)-A1)*CN1                                             WALL0048
      C1=ZWMAX-(RWMAX-A1)**2/B1                                         WALL0049
      DLZW=(ZWMAX-PW(60))/40.                                           WALL0050
      DO 20 I=61,137,2                                                  WALL0051
      PW(I+1)=PW(I-1)+DLZW                                              WALL0052
 20   PW(I)=A1+SQRT (B1*(PW(I+1)-C1))                                   WALL0053
      PW(139)=RWMAX                                                     WALL0054
      PW(140)=ZWMAX                                                     WALL0055
      J=140                                                             WALL0056
      GO TO 4                                                           WALL0057
C                                                                       WALL0058
C         ARC      (IWALL=3)                                            WALL0059
C                                                                       WALL0060
 3    RBR=((ZWMAX-PW(60))**2+(RWMAX-PW(59))**2)/                        WALL0061
     1    (2.*((ZWMAX-PW(60))*SN-(RWMAX-PW(59))*CS))                    WALL0062
      CN6=SN-(ZWMAX-PW(60))/RBR                                         WALL0063
      THER=ATAN (CN6/SQRT (1.-CN6**2))                                  WALL0064
      DLTHER=(THER-THJR)/40.                                            WALL0065
      ANG=THJR                                                          WALL0066
      DO 30 I=61,137,2                                                  WALL0067
      ANG=ANG+DLTHER                                                    WALL0068
      SNE=SIN (ANG)                                                     WALL0069
      CSE=COS (ANG)                                                     WALL0070
      PW(I)=PW(59)+RBR*(CSE-CS)                                         WALL0071
 30   PW(I+1)=PW(60)+RBR*(SN-SNE)                                       WALL0072
      PW(139)=RWMAX                                                     WALL0073
      PW(140)=ZWMAX                                                     WALL0074
      J=140                                                             WALL0075
C                                                                       WALL0076
 4    PW(J+1)=0.                                                        WALL0077
      PW(J+2)=0.                                                        WALL0078
      RETURN                                                            WALL0079
C                                                                       WALL0080
 999  WRITE (6,99)                                                      WALL0081
 99   FORMAT(19H0IWALL OUT OF RANGE )                                   WALL0082
      IERR=1                                                            WALL0083
      RETURN                                                            WALL0084
C                                                                       WALL0085
      END                                                               WALL0086
$ORIGIN        LEVEL2,UT2,REW
$IBFTC ACOMP                                                            ACOM0000
      SUBROUTINE ACOMP(M,N)                                             ACOM0001
      COMMON/NAMEP/P(50,100)                                            ACOM0002
      DIMENSION   A(10)    ,BS(10)   ,WPWG(10)           ,PP(50,100),   ACOM0003
     X PS(50,50) ,PTB(50,8),PW(150)  ,R(10)    ,TRL1(5)  ,TRL2(5)   ,   ACOM0004
     X TRL3(5)   ,TRL4(5)  ,TRM1(5)  ,TRM2(5)  ,TRM3(5)  ,TRM4(5)   ,   ACOM0005
     X FA(100)   ,GF(10)   ,FILL(24)                                    ACOM0006
      COMMON CAPN  ,CPL   ,CPS   ,CTOP  ,DL    ,DR    ,EK    ,EKW   ,   ACOM0007
     XEW    ,WPWG  ,GAMMA ,GMG0  ,HPL   ,HPS   ,IMAX  ,N1    ,N2    ,   ACOM0008
     XPS    ,PG0   ,PW    ,R     ,RCAP  ,RHOP0 ,SMP   ,TG0   ,TPM   ,   ACOM0009
     XA     ,BS    ,CC    ,CONST ,EOW1  ,EOW2  ,HIT   ,ICNT  ,IPP   ,   ACOM0010
     XJMAX  ,MW           ,PP    ,PTB   ,TRL1  ,TRL2  ,TRL3  ,TRL4  ,   ACOM0011
     XTRM1  ,TRM2  ,TRM3  ,TRM4  ,LRC   ,N3    ,FA    ,VGR   ,CPG   ,   ACOM0012
     XNILP  ,THI   ,GF    ,WDG   ,IFG   ,FILL  ,CG                      ACOM0013
C                                                                       ACOM0014
      IF(JMAX) 3,3,1                                                    ACOM0015
 1    DO 2 J=1,JMAX                                                     ACOM0016
      REJM=SQRT ((PTB(3,M)-PTB(J+30,M))**2 +(PTB(4,M)-PTB(J+40,M))**2)  ACOM0017
      REJN=SQRT ((PTB(3,N)-PTB(J+30,N))**2 +(PTB(4,N)-PTB(J+40,N))**2)  ACOM0018
      ARG=R(J)/CC*(REJM/PTB(7,M)*PTB(2,M)+REJN/PTB(7,N)*PTB(2,N))       ACOM0019
      CALL FIND(IFG,FA,ARG,FP)                                          ACOM0020
      CBAR=(REJM/SQRT (PTB(1,M)/PTB(2,M))+REJN/SQRT (PTB(1,N)/PTB(2,N)))ACOM0021
     X    /ARG                                                          ACOM0022
      GF(J)=FP/(1.+CG*CBAR*FP)                                          ACOM0023
      FF=FP*(((1.+9.45*CBAR)*(1.+2.52*CBAR)+3.03*CBAR**2)/              ACOM0024
     X ((1.+9.45*CBAR)*(1.+3.78*CBAR)+.573*(4.+1.134*CBAR)*CBAR**2))    ACOM0025
 2    A(J)=CONST/R(J)*FF/R(J)                                           ACOM0026
 3    RETURN                                                            ACOM0027
       END                                                              ACOM0028
$IBFTC ADJK                                                             ADJK0000
      SUBROUTINE ADJK(BARL,Q)                                           ADJK0001
C                                                                       ADJK0002
        COMMON/ERR/IERR                                                 ADJK0003
      COMMON/NAMEP/P(50,100)                                            ADJK0004
      DIMENSION   A(10)    ,BS(10)   ,WPWG(10)           ,PP(50,100),   ADJK0005
     X PS(50,50) ,PTB(50,8),PW(150)  ,R(10)    ,TRL1(5)  ,TRL2(5)   ,   ADJK0006
     X TRL3(5)   ,TRL4(5)  ,TRM1(5)  ,TRM2(5)  ,TRM3(5)  ,TRM4(5)  ,    ADJK0007
     X FA(100)   ,WDP(10)  ,GF(10)                                      ADJK0008
      COMMON CAPN  ,CPL   ,CPS   ,CTOP  ,DL    ,DR    ,INSERT,EKW   ,   ADJK0009
     XEW    ,WPWG  ,GAMMA ,GMG0  ,HPL   ,HPS   ,IMAX  ,N1    ,N2    ,   ADJK0010
     XPS    ,PG0   ,PW    ,R     ,RCAP  ,RHOP0 ,SMP   ,TG0   ,TPM   ,   ADJK0011
     XA     ,BS    ,CC    ,CONST ,EOW1  ,EOW2  ,HIT   ,ICNT  ,IPP   ,   ADJK0012
     XJMAX  ,MW           ,PP    ,PTB   ,TRL1  ,TRL2  ,TRL3  ,TRL4  ,   ADJK0013
     XTRM1  ,TRM2  ,TRM3  ,TRM4  ,LRC   ,N3    ,FA    ,VGR   ,CPG   ,   ADJK0014
     XNILP  ,THI   ,GF    ,WDG   ,IFG   ,RHOG0 ,THRUST,TI    ,WDP   ,   ADJK0015
     XTMD   ,IPPX                                                       ADJK0016
C                                                                       ADJK0017
      C67=0.                                                            ADJK0018
 1    C67=C67+1.                                                        ADJK0019
      IF(C67-7.) 3,85,85                                                ADJK0020
 3    DO 4 I=1,50                                                       ADJK0021
      PTB(I,6)=PP(I,IPP-1)                                              ADJK0022
 4    PTB(I,7)=PP(I,IPP)                                                ADJK0023
      STARK=(PTB(5,7)-PTB(5,6))/(PTB(6,7)-PTB(6,6))                     ADJK0024
      PTB(6,1)=(PTB(5,3)-PTB(5,6)+STARK*PTB(6,6)-BARL*PTB(6,3))         ADJK0025
     X         /(STARK-BARL)                                            ADJK0026
      Q=(PTB(6,1)-PTB(6,6))/(PTB(6,7)-PTB(6,6))                         ADJK0027
      IF(Q-1.) 5,5,7                                                    ADJK0028
 5    IF(Q) 6,999,999                                                   ADJK0029
 6    IPP=IPP-1                                                         ADJK0030
      IF(IPP-IPPX) 65,1,1                                               ADJK0031
 65   CALL ERROR(9)                                                     ADJK0032
      IERR=1                                                            ADJK0033
      RETURN                                                            ADJK0034
 7    IF(PP(1,IPP+1)) 8,9,8                                             ADJK0035
 8    IPP=IPP+1                                                         ADJK0036
      IF(PP(8,IPP-1)) 85,1,85                                           ADJK0037
 85   INSERT=1                                                          ADJK0038
      IPP=IPPX                                                          ADJK0039
      GO TO 999                                                         ADJK0040
 9    INSERT=-1                                                         ADJK0041
 999  RETURN                                                            ADJK0042
       END                                                              ADJK0043
$IBFTC AXISPT                                                           AXIS0000
      SUBROUTINE AXISPT                                                 AXIS0001
C     THE AXIS POINT COMPUTATION                                        AXIS0002
C                                                                       AXIS0003
      COMMON/NAMEP/P(50,100)                                            AXIS0004
      DIMENSION   A(10)    ,BS(10)   ,FJ(10)             ,PP(50,100),   AXIS0005
     X PS(50,50) ,PTB(50,8),PW(150)  ,R(10)    ,TRL1(5)  ,TRL2(5)   ,   AXIS0006
     X TRL3(5)   ,TRL4(5)  ,TRM1(5)  ,TRM2(5)  ,TRM3(5)  ,TRM4(5)  ,    AXIS0007
     X FA(100)                                                          AXIS0008
     X,RT1(3)    ,RT2(3)   ,RT3(3)   ,VAL1(3)  ,CAPL(3)                 AXIS0009
     X,Q(3)                                                             AXIS0010
      COMMON CAPN  ,CPL   ,CPS   ,CTOP  ,DL    ,DR    ,EK    ,EKW   ,   AXIS0011
     XEW    ,FJ    ,GAMMA ,GMG0  ,HPL   ,HPS   ,IMAX  ,N1    ,N2    ,   AXIS0012
     XPS    ,PG0   ,PW    ,R     ,RCAP  ,RHOP0 ,SMP   ,TG0   ,TPM   ,   AXIS0013
     XA     ,BS    ,CC    ,CONST ,EOW1  ,EOW2  ,HIT   ,ICNT  ,IPP   ,   AXIS0014
     XJMAX  ,MW           ,PP    ,PTB   ,TRL1  ,TRL2  ,TRL3  ,TRL4  ,   AXIS0015
     XTRM1  ,TRM2  ,TRM3  ,TRM4  ,LRC   ,DELTA ,FA    ,VGR              AXIS0016
C                                                                       AXIS0017
C     STEPS 1 AND 2 SETTING INITIAL CONDITIONS                          AXIS0018
C                                                                       AXIS0019
      PTB(1,3)=(PTB(1,1)+PTB(1,5))/2.                                   AXIS0020
      PTB(2,3)=(PTB(2,1)+PTB(2,5))/2.                                   AXIS0021
      PTB(3,3)=(PTB(3,1)+PTB(3,5))/2.                                   AXIS0022
      PTB(4,3)=0.                                                       AXIS0023
      PTB(5,3)=0.                                                       AXIS0024
      RT1(1)=SQRT (GAMMA*PTB(1,1))                                      AXIS0025
      RT2(1)=SQRT (PTB(2,1)*(PTB(3,1)**2 +PTB(4,1)**2 )-GAMMA*PTB(1,1)) AXIS0026
      RT3(1)=RT1(1)/RT2(1)                                              AXIS0027
      VAL1(1)=PTB(3,1)**2 *PTB(2,1)-GAMMA*PTB(1,1)                      AXIS0028
      GL1=(PTB(2,1)*PTB(3,1)*PTB(4,1)-RT1(1)*RT2(1))/VAL1(1)            AXIS0029
      GK1=(PTB(2,1)*PTB(3,1)*PTB(4,1)+RT1(1)*RT2(1))/VAL1(1)            AXIS0030
      ETA1=PTB(2,1)*(PTB(4,1)*RT1(1)-PTB(3,1)*RT2(1))/(RT1(1)*VAL1(1))  AXIS0031
      IF(JMAX) 2,2,11                                                   AXIS0032
 11   DO 19 J=1,JMAX                                                    AXIS0033
      PTB(J+10,3)=(PTB(J+10,1)+PTB(J+10,5))/2.                          AXIS0034
      PTB(J+20,3)=(PTB(J+20,1)+PTB(J+20,5))/2.                          AXIS0035
      PTB(J+30,3)=(PTB(J+30,1)+PTB(J+30,5))/2.                          AXIS0036
 19   PTB(J+40,3)=0.                                                    AXIS0037
      CALL ACOMP(1,5)                                                   AXIS0038
      CALL SUMP1(1)                                                     AXIS0039
      CALL SUMP1(5)                                                     AXIS0040
      F1  =RT3(1)*(PTB(7,1)*TRM4(1)-GAMMA*PTB(1,1)*PTB(4,1)/PTB(5,1))   AXIS0041
      S5=PTB(2,5)/(GAMMA*PTB(1,5))                                      AXIS0042
      T5=S5/PTB(3,5)*PTB(7,5)*TRM4(5)                                   AXIS0043
      GO TO 3                                                           AXIS0044
 2    F1  =-RT3(1)*GAMMA*PTB(1,1)*PTB(4,1)/PTB(5,1)                     AXIS0045
      G5=PTB(2,5)*PTB(3,5)                                              AXIS0046
      CAPL(1)=RT3(1)*PTB(2,1)*PTB(4,1)                                  AXIS0047
      S5=PTB(2,5)/(GAMMA*PTB(1,5))                                      AXIS0048
      T3=0.                                                             AXIS0049
      T5=0.                                                             AXIS0050
C                                                                       AXIS0051
C     STEP 3 CALCULATE BARL AND Z3                                      AXIS0052
C            BEGIN ITERATION LOOP                                       AXIS0053
C                                                                       AXIS0054
 3    DO 80 I=1,IMAX                                                    AXIS0055
      DO 210 N=1,4                                                      AXIS0056
 210  PTB(N,8)=PTB(N,3)                                                 AXIS0057
      RT1(3)=SQRT (GAMMA*PTB(1,3))                                      AXIS0058
      RT2(3)=SQRT (PTB(2,3)*(PTB(3,3)**2 +PTB(4,3)**2 )-GAMMA*PTB(1,3)) AXIS0059
      RT3(3)=RT1(3)/RT2(3)                                              AXIS0060
      VAL1(3)=PTB(3,3)**2 *PTB(2,3)-GAMMA*PTB(1,3)                      AXIS0061
      IF(JMAX) 23,23,211                                                AXIS0062
 211  DO 22 J=1,JMAX                                                    AXIS0063
      PTB(J+10,8)=PTB(J+10,3)                                           AXIS0064
      PTB(J+20,8)=PTB(J+20,3)                                           AXIS0065
      PTB(J+30,8)=PTB(J+30,3)                                           AXIS0066
 22   PTB(J+40,8)=PTB(J+40,3)                                           AXIS0067
 23   GL3=(PTB(2,3)*PTB(3,3)*PTB(4,3)-RT1(3)*RT2(3))/VAL1(3)            AXIS0068
      CALL TAFN(GL1,GL3,BARL)                                           AXIS0069
      PTB(6,3)=PTB(6,1)-PTB(5,1)/BARL                                   AXIS0070
      IF(JMAX)7,7,4                                                     AXIS0071
C                                                                       AXIS0072
C     STEP 4 CALCULATION OF PARTICLE H,U,AND RHO AT PT. 3               AXIS0073
C                                                                       AXIS0074
 4    PTB(7,3)=(PTB(1,3)/PTB(2,3))**CAPN                                AXIS0075
      C1=PTB(6,3)-PTB(6,5)                                              AXIS0076
      DO 49 J=1,JMAX                                                    AXIS0077
      B3=PTB(7,3)*(PTB(3,3)-PTB(J+30,3))/PTB(J+30,3)                    AXIS0078
      B5=PTB(7,5)*(PTB(3,5)-PTB(J+30,5))/PTB(J+30,5)                    AXIS0079
      PTB(J+30,3)=PTB(J+30,5)+A(J)*C1/2.*(B3+B5)                        AXIS0080
      PTB(J+10,3)=PTB(J+10,5)+C1/2.*(EFN(J,3)+EFN(J,5))                 AXIS0081
 49   PTB(J+20,3)=(PTB(J+20,5)*PTB(J+30,5)-2.*PTB(J+20,1)*PTB(J+40,1)   AXIS0082
     1   /PTB(5,1)*(PTB(6,3)-PTB(6,5)))/PTB(J+30,3)                     AXIS0083
C                                                                       AXIS0084
C     STEP 5 CALCULATION OF GAS PROPERTIES P, RHO, AND U                AXIS0085
C                                                                       AXIS0086
 7    ETA3=PTB(2,3)*(PTB(4,3)*RT1(3)-PTB(3,3)*RT2(3))/(RT1(3)*VAL1(3))  AXIS0087
      VGR =(PTB(4,1)/PTB(5,1))*PTB(3,3)/(PTB(3,1)+                      AXIS0088
     X     (PTB(6,3)-PTB(6,1))*PTB(4,1)/PTB(5,1))                       AXIS0089
      IF(JMAX) 75,75,71                                                 AXIS0090
 71   CALL SUMP1(3)                                                     AXIS0091
      F3  =RT3(3)*(PTB(7,3)*TRM4(3)-GAMMA*PTB(1,3)*VGR)                 AXIS0092
      C11=PTB(7,3)*TRM1(3)                                              AXIS0093
      C12=PTB(7,3)*TRM2(3)                                              AXIS0094
      C13=PTB(7,3)*TRM3(3)                                              AXIS0095
      CAPJ3=PTB(1,1)-(PTB(6,3)-PTB(6,1))*(F3  *ETA3+F1  *ETA1)          AXIS0096
     X      /2.-(RT3(1)+RT3(3))*(                                       AXIS0097
     X      PTB(4,1)*(PTB(2,1)*PTB(3,1)+PTB(2,3)*PTB(3,3))-             AXIS0098
     X      PTB(3,1)*(PTB(2,1)*PTB(4,1)+PTB(2,3)*PTB(4,3))              AXIS0099
     X      -(PTB(5,3)-PTB(5,1))*(C12+PTB(7,1)*TRM2(1))                 AXIS0100
     X      +(PTB(6,3)-PTB(6,1))*(C13+PTB(7,1)*TRM3(1))                 AXIS0101
     X      +(PTB(3,1)*(PTB(5,3)-PTB(5,1))-PTB(4,1)*(PTB(6,3)-PTB(6,1)) AXIS0102
     X      -PTB(3,3)*PTB(5,1)+PTB(4,3)*PTB(6,1))/2.*                   AXIS0103
     X      (C11+PTB(7,1)*TRM1(1)))/4.                                  AXIS0104
      CAPJ1=PTB(1,5)+((PTB(5,3)-PTB(5,5))*(C13+PTB(7,5)*TRM3(5))        AXIS0105
     X               +(PTB(6,3)-PTB(6,5))*(C12+PTB(7,5)*TRM2(5))        AXIS0106
     X   -PTB(3,5)/2.*(PTB(6,3)-PTB(6,5))*(C11+PTB(7,5)*TRM1(5))        AXIS0107
     X   -PTB(4,5)/2.*(PTB(5,3)-PTB(5,5))*(C11+PTB(7,5)*TRM1(5))        AXIS0108
     X      +PTB(4,5)*(PTB(2,3)*PTB(4,3)+PTB(2,5)*PTB(4,5))             AXIS0109
     X      +PTB(3,5)*(PTB(2,3)*PTB(3,3)+PTB(2,5)*PTB(3,5)))/2.         AXIS0110
      G3=PTB(2,3)*PTB(3,3)+(PTB(6,3)-PTB(6,5))*C11/2.                   AXIS0111
      G5=PTB(2,5)*PTB(3,5)+(PTB(6,3)-PTB(6,5))*PTB(7,5)*TRM1(5)/2.      AXIS0112
      CAPL(1)=RT3(1)*(PTB(2,1)*PTB(4,1)+PTB(5,3)*PTB(7,1)*TRM1(1)/2.)   AXIS0113
      S3=PTB(2,3)/(GAMMA*PTB(1,3))                                      AXIS0114
      T3=S3/PTB(3,3)*PTB(7,3)*TRM4(3)                                   AXIS0115
      GO TO 79                                                          AXIS0116
 75   F3  =-RT3(3)*GAMMA*PTB(1,3)*VGR                                   AXIS0117
      CAPJ3=PTB(1,1)-(PTB(6,3)-PTB(6,1))*(F3  *ETA3+F1  *ETA1)          AXIS0118
     X      /2.-(RT3(1)+RT3(3))*(                                       AXIS0119
     X      PTB(4,1)*(PTB(2,1)*PTB(3,1)+PTB(2,3)*PTB(3,3))-             AXIS0120
     X      PTB(3,1)*(PTB(2,1)*PTB(4,1)+PTB(2,3)*PTB(4,3)))/4.          AXIS0121
      CAPJ1=PTB(1,5)+(PTB(4,5)*(PTB(2,3)*PTB(4,3)+PTB(2,5)*             AXIS0122
     X    PTB(4,5))+PTB(3,5)*(PTB(2,3)*PTB(3,3)+PTB(2,5)*PTB(3,5)))/2.  AXIS0123
      G3=PTB(2,3)*PTB(3,3)                                              AXIS0124
      S3=PTB(2,3)/(GAMMA*PTB(1,3))                                      AXIS0125
 79   PTB(3,3)=2.*(CAPJ3-CAPJ1)/(CAPL(1)-G5-G3)                         AXIS0126
      PTB(1,3)=CAPJ3-CAPL(1)*PTB(3,3)/2.                                AXIS0127
      PTB(2,3)=PTB(2,5)+((S3+S5)*(PTB(1,3)-PTB(1,5))-(T3+T5)            AXIS0128
     X         *(PTB(6,3)-PTB(6,5)))/2.                                 AXIS0129
      CALL CHECK(DIFF)                                                  AXIS0130
      IF(DIFF) 80,81,80                                                 AXIS0131
 80   CONTINUE                                                          AXIS0132
      I=IMAX                                                            AXIS0133
 81   ICNT=I                                                            AXIS0134
      VGR =(PTB(4,1)/PTB(5,1))*PTB(3,3)/(PTB(3,1)+                      AXIS0135
     X     (PTB(6,3)-PTB(6,1))*PTB(4,1)/PTB(5,1))                       AXIS0136
      RETURN                                                            AXIS0137
       END                                                              AXIS0138
$IBFTC CHECK                                                            CHEC0000
      SUBROUTINE CHECK(DIFF)                                            CHEC0001
C                                                                       CHEC0002
      COMMON/NAMEP/P(50,100)                                            CHEC0003
      DIMENSION   A(10)    ,BS(10)   ,FJ(10)             ,PP(50,100),   CHEC0004
     X PS(50,50) ,PTB(50,8),PW(150)  ,R(10)    ,TRL1(5)  ,TRL2(5)   ,   CHEC0005
     X TRL3(5)   ,TRL4(5)  ,TRM1(5)  ,TRM2(5)  ,TRM3(5)  ,TRM4(5)       CHEC0006
      COMMON CAPN  ,CPL   ,CPS   ,CTOP  ,DL    ,DR    ,EK    ,EKW   ,   CHEC0007
     XEW    ,FJ    ,GAMMA ,GMG0  ,HPL   ,HPS   ,IMAX  ,N1    ,N2    ,   CHEC0008
     XPS    ,PG0   ,PW    ,R     ,RCAP  ,RHOP0 ,SMP   ,TG0   ,TPM   ,   CHEC0009
     XA     ,BS    ,CC    ,CONST ,EOW1  ,EOW2  ,HIT   ,ICNT  ,IPP   ,   CHEC0010
     XJMAX  ,MW           ,PP    ,PTB   ,TRL1  ,TRL2  ,TRL3  ,TRL4  ,   CHEC0011
     XTRM1  ,TRM2  ,TRM3  ,TRM4                                         CHEC0012
C                                                                       CHEC0013
      DO 1 I=1,4                                                        CHEC0014
      CALL CRIT(PTB(I,3),PTB(I,8),DIFF)                                 CHEC0015
      IF(DIFF) 5,1,5                                                    CHEC0016
 1    CONTINUE                                                          CHEC0017
      IF (JMAX) 5,5,3                                                   CHEC0018
 3    DO 2 I=1,4                                                        CHEC0019
      DO 2 J=1,JMAX                                                     CHEC0020
      JJ=10*I+J                                                         CHEC0021
      CALL CRIT(PTB(JJ,3),PTB(JJ,8),DIFF)                               CHEC0022
      IF(DIFF) 5,2,5                                                    CHEC0023
 2    CONTINUE                                                          CHEC0024
 5    RETURN                                                            CHEC0025
       END                                                              CHEC0026
$IBFTC CNTRL                                                            CNTR0000
      SUBROUTINE CNTRL                                                  CNTR0001
C                                                                       CNTR0002
        COMMON/ERR/IERR                                                 CNTR0003
      COMMON/NAMEP/P(50,100)                                            CNTR0004
      DIMENSION   A(10)    ,BS(10)   ,WPWG(10)           ,PP(50,100),   CNTR0005
     X PS(50,50) ,PTB(50,8),PW(150)  ,R(10)    ,TRL1(5)  ,TRL2(5)   ,   CNTR0006
     X TRL3(5)   ,TRL4(5)  ,TRM1(5)  ,TRM2(5)  ,TRM3(5)  ,TRM4(5)   ,   CNTR0007
     X FA(100)   ,WDP(10)  ,GF(10)                                      CNTR0008
      COMMON CAPN  ,CPL   ,CPS   ,CTOP  ,DL    ,DR    ,INSERT,DTWI  ,   CNTR0009
     XEW    ,WPWG  ,GAMMA ,GMG0  ,HPL   ,HPS   ,IMAX  ,N1    ,N2    ,   CNTR0010
     XPS    ,PG0   ,PW    ,R     ,RCAP  ,RHOP0 ,SMP   ,TG0   ,TPM   ,   CNTR0011
     XA     ,BS    ,CC    ,CONST ,EOW1  ,EOW2  ,HIT   ,ICNT  ,IPP   ,   CNTR0012
     XJMAX  ,MW           ,PP    ,PTB   ,TRL1  ,TRL2  ,TRL3  ,TRL4  ,   CNTR0013
     XTRM1  ,TRM2  ,TRM3  ,TRM4  ,LRC   ,DELTA ,FA    ,VGR   ,CPG   ,   CNTR0014
     XNILP  ,THI   ,GF    ,WDG   ,IFG   ,RHOG0 ,THRUST,TI    ,WDP   ,   CNTR0015
     XTMD   ,IPPX  ,BC    ,THII  ,DZMIN ,ITOT                           CNTR0016
     X,ZAX,RRT,RT,HOLD(116),IP,IT(20),KT(20),LAST                       CNTR0017
C                                                                       CNTR0018
C     AREA 0                                                            CNTR0019
C                                                                       CNTR0020
      STOP=0.                                                           CNTR0021
      DTW=DTWI*.017453293                                               CNTR0022
      CIMAX=15.                                                         CNTR0023
      CLMAX=15.                                                         CNTR0024
      CRMAX=15.                                                         CNTR0025
      CWMAX=15.                                                         CNTR0026
      DL2=DL*DL                                                         CNTR0027
      DR2=DR*DR                                                         CNTR0028
      IPMX=99                                                           CNTR0029
      LASTMX=20                                                         CNTR0030
C                                                                       CNTR0031
C     AREA 1                                                            CNTR0032
C                                                                       CNTR0033
      JMAXI=JMAX                                                        CNTR0034
      MW=4                                                              CNTR0035
      MS=1                                                              CNTR0036
      LAST=1                                                            CNTR0037
      IT(1)=0                                                           CNTR0038
      EOW1=0.                                                           CNTR0039
      EOW2=0.                                                           CNTR0040
      INSERT=0                                                          CNTR0041
      DO 1 I=1,50                                                       CNTR0042
      PTB(I,1)=PS(I,1)                                                  CNTR0043
 1    PP(I,1)=PS(I,1)                                                   CNTR0044
      LRC=0                                                             CNTR0045
      CALL PRINT (1)                                                    CNTR0046
      LRC=1                                                             CNTR0047
 2    MS=MS+1                                                           CNTR0048
      IP=1                                                              CNTR0049
      IPP=1                                                             CNTR0050
      DO 3 I=1,50                                                       CNTR0051
      P(I,1)=PS(I,MS)                                                   CNTR0052
      PTB(I,1)=PP(I,1)                                                  CNTR0053
 3    PTB(I,2)=P(I,1)                                                   CNTR0054
      CALL PRINT(1)                                                     CNTR0055
      CI=0.                                                             CNTR0056
      ID=2                                                              CNTR0057
 4    CALL PTINT(0)                                                     CNTR0058
      IF((PTB(5,3)-PTB(5,2))**2+(PTB(6,3)-PTB(6,2))**2-DL2) 5,5,8       CNTR0059
 5    IF((PTB(5,3)-PTB(5,1))**2+(PTB(6,3)-PTB(6,1))**2-DR2) 6,6,8       CNTR0060
 6    IF(CI) 41,12,41                                                   CNTR0061
 8    CI=CI+1.                                                          CNTR0062
      IF(CI-CIMAX) 10,9,9                                               CNTR0063
 9    CALL ERROR(1)                                                     CNTR0064
      RETURN                                                            CNTR0065
 10   DO 11 I=1,50                                                      CNTR0066
 11   PTB(I,2)=(PP(I,1)+PTB(I,2))/2.                                    CNTR0067
      PTB(7,2)=(PTB(1,2)/PTB(2,2))**CAPN                                CNTR0068
      GO TO 4                                                           CNTR0069
C                                                                       CNTR0070
C     AREA 2                                                            CNTR0071
C                                                                       CNTR0072
 12   IPP=IPP+1                                                         CNTR0073
 13   IP=IP+1                                                           CNTR0074
      CALL PRINT(ID)                                                    CNTR0075
      IF(STOP) 130,132,130                                              CNTR0076
 130  IF(IP-IPMX) 137,131,131                                           CNTR0077
 131  CALL ERROR(12)                                                    CNTR0078
      RETURN                                                            CNTR0079
 132  IF(IP-IPMX+10) 137,137,133                                        CNTR0080
 133  DL2=1.                                                            CNTR0081
      DR2=1.                                                            CNTR0082
      DTW=1.                                                            CNTR0083
      STOP=1.                                                           CNTR0084
 137  IF(PTB(8,3)) 138,139,138                                          CNTR0085
 138  JMAX=JMAX-1                                                       CNTR0086
 139  DO 14 I=1,50                                                      CNTR0087
      P(I,IP)=PTB(I,3)                                                  CNTR0088
 14   PTB(I,3)=0.                                                       CNTR0089
      IF(EOW2) 999,15,999                                               CNTR0090
C                                                                       CNTR0091
C     AREA 3                                                            CNTR0092
C                                                                       CNTR0093
 15   IF(PP(1,IPP-1)) 35,16,35                                          CNTR0094
 16   IF(IT(LAST)) 17,21,17                                             CNTR0095
 17   PP(1,IP+1)=0.                                                     CNTR0096
 171  II=IT(LAST)                                                       CNTR0097
      DO 18 J=II,IP                                                     CNTR0098
      DO 18 I=1,50                                                      CNTR0099
 18   PP(I,J)=P(I,J)                                                    CNTR0100
      IP=IT(LAST)-1                                                     CNTR0101
      IPP=IT(LAST)+1                                                    CNTR0102
      JMAX=KT(LAST)                                                     CNTR0103
      LAST=LAST-1                                                       CNTR0104
      GO TO 49                                                          CNTR0105
C                                                                       CNTR0106
C     AREA 4                                                            CNTR0107
C                                                                       CNTR0108
 21   PP(1,IP+1)=0.                                                     CNTR0109
 210  DO 22 J=1,IP                                                      CNTR0110
      DO 22 I=1,50                                                      CNTR0111
      PP(I,J)=P(I,J)                                                    CNTR0112
 22   P(I,J)=0.                                                         CNTR0113
      CK=0.                                                             CNTR0114
      LRC=LRC+1                                                         CNTR0115
      IF(PS(1,MS+1)) 23,26,23                                           CNTR0116
 23   IF(PS(8,MS)) 24,25,24                                             CNTR0117
 24   JMAXI=JMAXI+1                                                     CNTR0118
 25   JMAX=JMAXI                                                        CNTR0119
      GO TO 2                                                           CNTR0120
C                                                                       CNTR0121
C     AREA 5                                                            CNTR0122
C                                                                       CNTR0123
 26   IP=0                                                              CNTR0124
      IPP=2                                                             CNTR0125
      JMAX=JMAXI                                                        CNTR0126
      DO 27 I=1,50                                                      CNTR0127
      PTB(I,5)=PP(I,1)                                                  CNTR0128
 27   PTB(I,1)=PP(I,2)                                                  CNTR0129
      CR=0.                                                             CNTR0130
      ID=3                                                              CNTR0131
 271  CALL AXISPT                                                       CNTR0132
      IF((PTB(5,3)-PTB(5,1))**2+(PTB(6,3)-PTB(6,1))**2-DR2) 28,28,31    CNTR0133
 28   IF(CR) 30,12,30                                                   CNTR0134
 30   CALL PRINT(6)                                                     CNTR0135
      GO TO 13                                                          CNTR0136
 31   CR=CR+1.                                                          CNTR0137
      IF(CR-CRMAX) 33,32,32                                             CNTR0138
 32   CALL ERROR(2)                                                     CNTR0139
      RETURN                                                            CNTR0140
 33   DO 34 I=1,50                                                      CNTR0141
 34   PTB(I,1)=(PTB(I,5)+PTB(I,1))/2.                                   CNTR0142
      PTB(7,1)=(PTB(1,1)/PTB(2,1))**CAPN                                CNTR0143
      GO TO 271                                                         CNTR0144
C                                                                       CNTR0145
C     AREA 6                                                            CNTR0146
C                                                                       CNTR0147
 35   IF(PP(1,IPP)) 49,37,49                                            CNTR0148
 37   DO 38 I=1,50                                                      CNTR0149
      PTB(I,5)=PP(I,IPP-1)                                              CNTR0150
 38   PTB(I,2)=P(I,IP)                                                  CNTR0151
      CR=0.                                                             CNTR0152
      ID=5                                                              CNTR0153
 381  CALL WLPT                                                         CNTR0154
      IF(IERR .NE. 0) RETURN                                            CNTR0155
      IF(EOW1) 45,39,45                                                 CNTR0156
 39   IF(ABS (ATAN (PTB(4,5)/PTB(3,5))-ATAN (PTB(4,3)/PTB(3,3)))        CNTR0157
     1  -DTW) 40,40,45                                                  CNTR0158
 40   T1=PTB(5,3)-PTB(5,5)                                              CNTR0159
      T2=.04882191*(PTB(5,3)+PTB(5,5))                                  CNTR0160
      DT=(PTB(1,3)+PTB(1,5))*T2*T1                                      CNTR0161
      THRUST=THRUST+DT*RT*RT                                            CNTR0162
 401  IF(CR) 41,12 ,41                                                  CNTR0163
 41   CALL PRINT(7)                                                     CNTR0164
      IP=IP+1                                                           CNTR0165
      DO 42 I=1,50                                                      CNTR0166
 42   P(I,IP)=PTB(I,2)                                                  CNTR0167
      P(8,IP)=0.                                                        CNTR0168
      LAST=LAST+1                                                       CNTR0169
      IF(LAST-LASTMX) 44,44,43                                          CNTR0170
 43   CALL ERROR(6)                                                     CNTR0171
      RETURN                                                            CNTR0172
 44   IT(LAST)=IP                                                       CNTR0173
      KT(LAST)=JMAX                                                     CNTR0174
      GO TO 12                                                          CNTR0175
 45   CR=CR+1.                                                          CNTR0176
      IF(CR-CWMAX) 47,46,46                                             CNTR0177
 46   CALL ERROR(3)                                                     CNTR0178
      RETURN                                                            CNTR0179
 47   DO 48 I=1,50                                                      CNTR0180
 48   PTB(I,2)=(PTB(I,5)+PTB(I,2))/2.                                   CNTR0181
      PTB(7,2)=(PTB(1,2)/PTB(2,2))**CAPN                                CNTR0182
      EOW2=0.                                                           CNTR0183
      GO TO 381                                                         CNTR0184
C                                                                       CNTR0185
C     AREA 7                                                            CNTR0186
C                                                                       CNTR0187
 49   CR=0.                                                             CNTR0188
      IPPX=IPP                                                          CNTR0189
      IF (PP(8,IPP-1)) 50,52,50                                         CNTR0190
 50   DO 51 I=1,50                                                      CNTR0191
      PTB(I,4)=PP(I,IPP-1)                                              CNTR0192
 51   PTB(I,2)=P(I,IP)                                                  CNTR0193
      CL=-1.                                                            CNTR0194
      ID=4                                                              CNTR0195
      ASSIGN 511 TO N                                                   CNTR0196
 511  INSERT=0                                                          CNTR0197
      CALL KPBPT                                                        CNTR0198
      IF(IERR .NE. 0) RETURN                                            CNTR0199
      IPP=IPP                                                           CNTR0200
      IF(INSERT) 571,57,571                                             CNTR0201
C                                                                       CNTR0202
 52   IPPP=IPP-1                                                        CNTR0203
 521  IF (PP(5,IPP)-PP(5,IPPP)-DELTA) 522,529,529                       CNTR0204
 522  IF(PP(8,IPP)) 529,523,529                                         CNTR0205
 523  IF(PP(1,IPP+1)) 524,529,524                                       CNTR0206
 524  IPP=IPP+1                                                         CNTR0207
      GO TO 521                                                         CNTR0208
 529  DO 53 I=1,50                                                      CNTR0209
      PTB(I,4)=PP(I,IPPP)                                               CNTR0210
      PTB(I,1)=PP(I,IPP)                                                CNTR0211
 53   PTB(I,2)=P(I,IP)                                                  CNTR0212
      CL=0.                                                             CNTR0213
      ID=2                                                              CNTR0214
      ASSIGN 531 TO N                                                   CNTR0215
 531  CALL PTINT(4)                                                     CNTR0216
      IF(PTB(6,3)-PP(6,IPP)+(PP(5,IPP)-PTB(5,3))/(PP(5,IPP)-            CNTR0217
     X  PP(5,IPP-1))*(PP(6,IPP)-PP(6,IPP-1))) 70,70,533                 CNTR0218
 533  IF((PTB(5,3)-PTB(5,2))**2+(PTB(6,3)-PTB(6,2))**2-DL2) 57,57,54    CNTR0219
 54   CL=CL+1.                                                          CNTR0220
      IF(CL-CLMAX) 55,541,541                                           CNTR0221
 541  CALL ERROR(5)                                                     CNTR0222
      RETURN                                                            CNTR0223
 55   DO 56 I=1,50                                                      CNTR0224
 56   PTB(I,1)=(PP(I,IPPX-1)+PTB(I,1))/2.                               CNTR0225
      PTB(7,1)=(PTB(1,1)/PTB(2,1))**CAPN                                CNTR0226
      GO TO 531                                                         CNTR0227
 57   IF((PTB(5,3)-PTB(5,1))**2+(PTB(6,3)-PTB(6,1))**2-DR2) 61,61,571   CNTR0228
 571  CR=CR+1.                                                          CNTR0229
      IF(CR-CRMAX) 59,58,58                                             CNTR0230
 58   CALL ERROR(4)                                                     CNTR0231
      RETURN                                                            CNTR0232
 59   DO 60 I=1,50                                                      CNTR0233
 60   PTB(I,2)=(PP(I,IPPX-1)+PTB(I,2))/2.                               CNTR0234
      PTB(7,2)=(PTB(1,2)/PTB(2,2))**CAPN                                CNTR0235
      GO TO N, (511,531)                                                CNTR0236
 61   IF(CL) 62,401,62                                                  CNTR0237
 62   IPP=IPP-1                                                         CNTR0238
      CALL PRINT(6)                                                     CNTR0239
      DO 63 I=1,50                                                      CNTR0240
 63   PP(I,IPP)=PTB(I,1)                                                CNTR0241
      PP(8,IPP)=0.                                                      CNTR0242
      GO TO 401                                                         CNTR0243
C                                                                       CNTR0244
C     AREA 8                                                            CNTR0245
C                                                                       CNTR0246
 70   WRITE (6,700)                                                     CNTR0247
 700  FORMAT(6H SHOCK)                                                  CNTR0248
      KK=IPP-IP-1                                                       CNTR0249
      IF(KK) 710,76,73                                                  CNTR0250
 710  DO 711 JJ=IPP,100                                                 CNTR0251
      IF(PP(1,JJ)) 71,72,711                                            CNTR0252
 711  CONTINUE                                                          CNTR0253
 71   WRITE (6,701)                                                     CNTR0254
 701  FORMAT(12H SHOCK ERROR)                                           CNTR0255
      IERR=1                                                            CNTR0256
      RETURN                                                            CNTR0257
 72   JJ=JJ                                                             CNTR0258
      JPP=JJ+IPP                                                        CNTR0259
      DO 720 J=IPP,JJ                                                   CNTR0260
      I=JPP-J                                                           CNTR0261
      K=I-KK                                                            CNTR0262
      DO 720 L=1,50                                                     CNTR0263
 720  PP(L,K)=PP(L,I)                                                   CNTR0264
      GO TO 76                                                          CNTR0265
 73   DO 75 J=IPP,100                                                   CNTR0266
      K=J-KK                                                            CNTR0267
      DO 74 I=1,50                                                      CNTR0268
 74   PP(I,K)=PP(I,J)                                                   CNTR0269
      IF(PP(1,J)) 71,76,75                                              CNTR0270
 75   CONTINUE                                                          CNTR0271
      GO TO 71                                                          CNTR0272
 76   IF(IT(LAST)) 71,210,171                                           CNTR0273
 999  RETURN                                                            CNTR0274
       END                                                              CNTR0275
$IBFTC CRIT                                                             CRIT0000
      SUBROUTINE CRIT(XN,XM,IDIF)                                       CRIT0001
      IDIF=0                                                            CRIT0002
      IF(XN .EQ. 0.) RETURN                                             CRIT0003
      IF(ABS((XM-XN)/XN).GT.5.0E-5) IDIF=1                              CRIT0004
      RETURN                                                            CRIT0005
      END                                                               CRIT0006
$IBFTC EFN                                                              EFN 0000
      FUNCTION EFN(J,N)                                                 EFN 0001
C                                                                       EFN 0002
      COMMON/NAMEP/P(50,100)                                            EFN 0003
      DIMENSION   A(10)    ,BS(10)   ,WPWG(10)           ,PP(50,100),   EFN 0004
     X PS(50,50) ,PTB(50,8),PW(150)  ,R(10)    ,TRL1(5)  ,TRL2(5)   ,   EFN 0005
     X TRL3(5)   ,TRL4(5)  ,TRM1(5)  ,TRM2(5)  ,TRM3(5)  ,TRM4(5)   ,   EFN 0006
     X FA(100)   ,GF(10)                                                EFN 0007
      COMMON CAPN  ,CPL   ,CPS   ,CTOP  ,DL    ,DR    ,EK    ,EKW   ,   EFN 0008
     XEW    ,WPWG  ,GAMMA ,GMG0  ,HPL   ,HPS   ,IMAX  ,N1    ,N2    ,   EFN 0009
     XPS    ,PG0   ,PW    ,R     ,RCAP  ,RHOP0 ,SMP   ,TG0   ,TPM   ,   EFN 0010
     XA     ,BS    ,CC    ,CONST ,EOW1  ,EOW2  ,HIT   ,ICNT  ,IPP   ,   EFN 0011
     XJMAX  ,MW           ,PP    ,PTB   ,TRL1  ,TRL2  ,TRL3  ,TRL4  ,   EFN 0012
     XTRM1  ,TRM2  ,TRM3  ,TRM4  ,LRC   ,N3    ,FA    ,VGR   ,CPG   ,   EFN 0013
     XNILP  ,THI   ,GF                                                  EFN 0014
C                                                                       EFN 0015
      C1=BS(J)*PTB(7,N)/PTB(J+30,N)*GF(J)                               EFN 0016
      C2=PTB(1,N) /RCAP/PTB(2,N)                                        EFN 0017
      IF (PTB(J+10,N)-HPL) 10,10,30                                     EFN 0018
 10   IF (PTB(J+10,N)-HPS) 11,20,20                                     EFN 0019
 11   EFN=-C1*(PTB(J+10,N)/CPS-C2)                                      EFN 0020
      GO TO 40                                                          EFN 0021
 20   EFN=-C1*(TPM-C2)                                                  EFN 0022
      GO TO 40                                                          EFN 0023
 30   EFN=-C1*(TPM+(PTB(J+10,N)-HPL)/CPL-C2)                            EFN 0024
 40   RETURN                                                            EFN 0025
       END                                                              EFN 0026
$IBFTC ERROR                                                            ERRO0000
      SUBROUTINE ERROR (I)                                              ERRO0001
C                                                                       ERRO0002
        COMMON/ERR/IERR                                                 ERRO0003
      COMMON/NAMEP/P(50,100)                                            ERRO0004
C                                                                       ERRO0005
      GO TO (1,2,3,4,5,6,7,8,9,10,11,12),I                              ERRO0006
 1    WRITE (6,101)                                                     ERRO0007
 101  FORMAT (48H0 ERROR(1). UNABLE TO INSERT ALONG INITIAL LINE.)      ERRO0008
      GO TO 999                                                         ERRO0009
 2    WRITE (6,102)                                                     ERRO0010
 102  FORMAT (44H0 ERROR(2). UNABLE TO INSERT FOR AXIS POINT.)          ERRO0011
      GO TO 999                                                         ERRO0012
 3    WRITE (6,103)                                                     ERRO0013
 103  FORMAT (44H0 ERROR(3). UNABLE TO INSERT FOR WALL POINT.)          ERRO0014
      GO TO 999                                                         ERRO0015
 4    WRITE (6,104)                                                     ERRO0016
 104  FORMAT (29H0 ERROR(4). UNABLE TO INSERT.)                         ERRO0017
      GO TO 999                                                         ERRO0018
 5    WRITE (6,105)                                                     ERRO0019
 105  FORMAT (29H0 ERROR(5). UNABLE TO INSERT.)                         ERRO0020
      GO TO 999                                                         ERRO0021
 6    WRITE (6,106)                                                     ERRO0022
 106  FORMAT (42H0 ERROR(6). TOO MANY INSERTIONS ATTEMPTED.)            ERRO0023
      GO TO 999                                                         ERRO0024
 7    WRITE (6,107)                                                     ERRO0025
 107  FORMAT (49H0 ERROR(7). KTH PART BNDRY PT UNABLE TO HIT WALL.)     ERRO0026
      GO TO 999                                                         ERRO0027
 8    WRITE (6,108)                                                     ERRO0028
 108  FORMAT (54H0 ERROR(8). UNABLE TO FIND PTS 6+7 FOR KTH PART BNDRY.)ERRO0029
      GO TO 999                                                         ERRO0030
 9    WRITE (6,109)                                                     ERRO0031
 109  FORMAT (43H0 ERROR(9). KTH PART BNDRY PT 1 BELOW PT 4.)           ERRO0032
      GO TO 999                                                         ERRO0033
 10   WRITE (6,110)                                                     ERRO0034
 110  FORMAT (46H0 ERROR(10) UNABLE TO FIND PTS 6+7 ALONG WALL.)        ERRO0035
      GO TO 999                                                         ERRO0036
 11   WRITE (6,111)                                                     ERRO0037
 111  FORMAT (42H0 ERROR(11) CHARACTORISTICS CROSS AT WALL.)            ERRO0038
      GO TO 999                                                         ERRO0039
 12   WRITE (6,112)                                                     ERRO0040
 112  FORMAT (43H0 ERROR(12) TOO MANY CHARACTORISTIC POINTS.)           ERRO0041
 999  IERR=1                                                            ERRO0042
      RETURN                                                            ERRO0043
       END                                                              ERRO0044
$IBFTC KPBPT                                                            KPBP0000
      SUBROUTINE KPBPT                                                  KPBP0001
C     THE KTH PARTICLE BOUNDARY CALCULATION                             KPBP0002
C                                                                       KPBP0003
        COMMON/ERR/IERR                                                 KPBP0004
      COMMON/NAMEP/P(50,100)                                            KPBP0005
      DIMENSION   A(10)    ,BS(10)   ,WPWG(10)           ,PP(50,100),   KPBP0006
     X PS(50,50) ,PTB(50,8),PW(150)  ,R(10)    ,TRL1(5)  ,TRL2(5)   ,   KPBP0007
     X TRL3(5)   ,TRL4(5)  ,TRM1(5)  ,TRM2(5)  ,TRM3(5)  ,TRM4(5)   ,   KPBP0008
     X FA(100)                                                          KPBP0009
     X,RT1(3)    ,RT2(3)   ,RT3(3)   ,VAL1(3)  ,CAPL(3)                 KPBP0010
     X,Q(3)      ,SAVE(7)                                               KPBP0011
     X,RHO4(10),U4(10),V4(10)                                           KPBP0012
      COMMON CAPN  ,CPL   ,CPS   ,CTOP  ,DL    ,DR    ,INSERT,EKW   ,   KPBP0013
     XEW    ,WPWG  ,GAMMA ,GMG0  ,HPL   ,HPS   ,IMAX  ,N1    ,N2    ,   KPBP0014
     XPS    ,PG0   ,PW    ,R     ,RCAP  ,RHOP0 ,SMP   ,TG0   ,TPM   ,   KPBP0015
     XA     ,BS    ,CC    ,CONST ,EOW1  ,EOW2  ,HIT   ,ICNT  ,IPP   ,   KPBP0016
     XJMAX  ,MW           ,PP    ,PTB   ,TRL1  ,TRL2  ,TRL3  ,TRL4  ,   KPBP0017
     XTRM1  ,TRM2  ,TRM3  ,TRM4  ,LRC   ,DELTA ,FA    ,VGR              KPBP0018
C                                                                       KPBP0019
C     STEP 1 AVERAGEING FOR INITIAL CONDITIONS                          KPBP0020
C                                                                       KPBP0021
      DO 1 I=1,4                                                        KPBP0022
      PTB(I,3)=(PTB(I,4)+PTB(I,2))/2.                                   KPBP0023
 1    PTB(I,1)=PTB(I,4)                                                 KPBP0024
      PTB(JMAX+30,1)=PTB(JMAX+30,2)                                     KPBP0025
      PTB(JMAX+40,1)=PTB(JMAX+40,2)                                     KPBP0026
      PTB(3,5)=PTB(3,3)                                                 KPBP0027
      PTB(4,5)=PTB(4,3)                                                 KPBP0028
      DO 12 N=1,7                                                       KPBP0029
 12   SAVE(N)=PTB(N,4)                                                  KPBP0030
      L=JMAX-1                                                          KPBP0031
      R4=PTB(5,4)                                                       KPBP0032
      Z4=PTB(6,4)                                                       KPBP0033
      DO 13 J=1,JMAX                                                    KPBP0034
      RHO4(J)=PTB(J+20,4)                                               KPBP0035
      U4(J)=PTB(J+30,4)                                                 KPBP0036
      V4(J)=PTB(J+40,4)                                                 KPBP0037
      PTB(J+10,3)=(PTB(J+10,4)+PTB(J+10,2))/2.                          KPBP0038
      PTB(J+20,3)=(PTB(J+20,4)+PTB(J+20,2))/2.                          KPBP0039
      PTB(J+30,3)=(PTB(J+30,4)+PTB(J+30,2))/2.                          KPBP0040
 13   PTB(J+40,3)=(PTB(J+40,4)+PTB(J+40,2))/2.                          KPBP0041
      RT1(2)=SQRT (GAMMA*PTB(1,2))                                      KPBP0042
      RT2(2)=SQRT (PTB(2,2)*(PTB(3,2)**2 +PTB(4,2)**2 )-GAMMA*PTB(1,2)) KPBP0043
      RT3(2)=RT1(2)/RT2(2)                                              KPBP0044
      VAL1(2)=PTB(3,2)**2 *PTB(2,2)-GAMMA*PTB(1,2)                      KPBP0045
      GK2=(PTB(2,2)*PTB(3,2)*PTB(4,2)+RT1(2)*RT2(2))/VAL1(2)            KPBP0046
      GMU2=PTB(2,2)*(PTB(4,2)*RT1(2)+PTB(3,2)*RT2(2))/(RT1(2)*VAL1(2))  KPBP0047
      IF (L) 2,2,14                                                     KPBP0048
 14   DO 16 J=1,L                                                       KPBP0049
      PTB(J+30,4)=PTB(J+30,3)                                           KPBP0050
 16   PTB(J+40,4)=PTB(J+40,3)                                           KPBP0051
C                                                                       KPBP0052
C     STEP 2 CALCULATION OF BARK, BARC, Z3,AND R3                       KPBP0053
C            BEGIN ITERATION LOOP                                       KPBP0054
C                                                                       KPBP0055
 2    DO 80 I=1,IMAX                                                    KPBP0056
      DO 21 N=1,3,2                                                     KPBP0057
      RT1(N)=SQRT (GAMMA*PTB(1,N))                                      KPBP0058
      RT2(N)=SQRT (PTB(2,N)*(PTB(3,N)**2 +PTB(4,N)**2 )-GAMMA*PTB(1,N)) KPBP0059
      RT3(N)=RT1(N)/RT2(N)                                              KPBP0060
 21   VAL1(N)=PTB(3,N)**2 *PTB(2,N)-GAMMA*PTB(1,N)                      KPBP0061
      DO 210 N=1,4                                                      KPBP0062
 210  PTB(N,8)=PTB(N,3)                                                 KPBP0063
      DO 22 J=1,JMAX                                                    KPBP0064
      PTB(J+10,8)=PTB(J+10,3)                                           KPBP0065
      PTB(J+20,8)=PTB(J+20,3)                                           KPBP0066
      PTB(J+30,8)=PTB(J+30,3)                                           KPBP0067
 22   PTB(J+40,8)=PTB(J+40,3)                                           KPBP0068
      GK3=(PTB(2,3)*PTB(3,3)*PTB(4,3)+RT1(3)*RT2(3))/VAL1(3)            KPBP0069
      CALL TAFN(GK2,GK3,BARK)                                           KPBP0070
      CALL TAFN(PTB(JMAX+40,3)/PTB(JMAX+30,3),PTB(JMAX+40,4)/           KPBP0071
     XPTB(JMAX+30,4),BARC)                                              KPBP0072
      PTB(6,3)=(PTB(5,4)-PTB(5,2)+BARK*PTB(6,2)-BARC*PTB(6,4))          KPBP0073
     X         /(BARK-BARC)                                             KPBP0074
      PTB(5,3)=PTB(5,2)+BARK*(PTB(6,3)-PTB(6,2))                        KPBP0075
      PTB(7,3)=(PTB(1,3)/PTB(2,3))**CAPN                                KPBP0076
C                                                                       KPBP0077
C     STEP 3 CALCULATE BARL                                             KPBP0078
C                                                                       KPBP0079
 31   GL1=(PTB(2,1)*PTB(3,1)*PTB(4,1)-RT1(1)*RT2(1))/VAL1(1)            KPBP0080
      GL3=(PTB(2,3)*PTB(3,3)*PTB(4,3)-RT1(3)*RT2(3))/VAL1(3)            KPBP0081
      CALL TAFN(GL1,GL3,BARL)                                           KPBP0082
C                                                                       KPBP0083
C     STEP 4 CALCULATE STARK, Z1, Q, AND DECIDE PTS. 6 AND 7            KPBP0084
C                                                                       KPBP0085
      CALL ADJK(BARL,Q)                                                 KPBP0086
      IF(IERR .NE. 0) RETURN                                            KPBP0087
      IF(INSERT) 81,5,81                                                KPBP0088
C                                                                       KPBP0089
C     STEP 5 CALCULATION OF THE POINT 1 BUCKET                          KPBP0090
C                                                                       KPBP0091
 5    PTB(5,1)=PTB(5,6)+Q*(PTB(5,7)-PTB(5,6))                           KPBP0092
      PTB(1,1)=PTB(1,6)+Q*(PTB(1,7)-PTB(1,6))                           KPBP0093
      PTB(2,1)=PTB(2,6)+Q*(PTB(2,7)-PTB(2,6))                           KPBP0094
      PTB(3,1)=PTB(3,6)+Q*(PTB(3,7)-PTB(3,6))                           KPBP0095
      PTB(4,1)=PTB(4,6)+Q*(PTB(4,7)-PTB(4,6))                           KPBP0096
      IF(L) 59,59,51                                                    KPBP0097
 51   DO 58 J=1,L                                                       KPBP0098
      PTB(J+10,1)=PTB(J+10,6)+Q*(PTB(J+10,7)-PTB(J+10,6))               KPBP0099
      PTB(J+20,1)=PTB(J+20,6)+Q*(PTB(J+20,7)-PTB(J+20,6))               KPBP0100
      PTB(J+30,1)=PTB(J+30,6)+Q*(PTB(J+30,7)-PTB(J+30,6))               KPBP0101
 58   PTB(J+40,1)=PTB(J+40,6)+Q*(PTB(J+40,7)-PTB(J+40,6))               KPBP0102
 59   PTB(7,1)=(PTB(1,1)/PTB(2,1))**CAPN                                KPBP0103
C                                                                       KPBP0104
C     STEPS 6,7, AND 8 CALCULATION OF THE J VARABLES                    KPBP0105
C                                                                       KPBP0106
 6    CALL ACOMP(1,2)                                                   KPBP0107
      DO 89 J=1,JMAX                                                    KPBP0108
      IF(J-JMAX) 61,70,70                                               KPBP0109
 61   CALL TAFN(PTB(J+40,3)/PTB(J+30,3),PTB(J+40,4)/PTB(J+30,4),BARC)   KPBP0110
      PTB(5,4)=(PTB(5,3)+BARC*(PTB(6,2)-PTB(6,3)-PTB(5,2)               KPBP0111
     X         *(PTB(6,1)-PTB(6,2))/(PTB(5,1)-PTB(5,2))))               KPBP0112
     X         /(1.-BARC*(PTB(6,1)-PTB(6,2))/(PTB(5,1)-PTB(5,2)))       KPBP0113
      Q=(PTB(5,4)-PTB(5,2))/(PTB(5,1)-PTB(5,2))                         KPBP0114
      PTB(6,4)=PTB(6,2)+Q*(PTB(6,1)-PTB(6,2))                           KPBP0115
      PTB(1,4)=PTB(1,2)+Q*(PTB(1,1)-PTB(1,2))                           KPBP0116
      PTB(2,4)=PTB(2,2)+Q*(PTB(2,1)-PTB(2,2))                           KPBP0117
      PTB(3,4)=PTB(3,2)+Q*(PTB(3,1)-PTB(3,2))                           KPBP0118
      PTB(4,4)=PTB(4,2)+Q*(PTB(4,1)-PTB(4,2))                           KPBP0119
      PTB(7,4)=(PTB(1,4)/PTB(2,4))**CAPN                                KPBP0120
      PTB(J+10,4)=PTB(J+10,2)+Q*(PTB(J+10,1)-PTB(J+10,2))               KPBP0121
      PTB(J+20,4)=PTB(J+20,2)+Q*(PTB(J+20,1)-PTB(J+20,2))               KPBP0122
      PTB(J+30,4)=PTB(J+30,2)+Q*(PTB(J+30,1)-PTB(J+30,2))               KPBP0123
      PTB(J+40,4)=PTB(J+40,2)+Q*(PTB(J+40,1)-PTB(J+40,2))               KPBP0124
      GO TO 72                                                          KPBP0125
 70   DO 71 N=1,7                                                       KPBP0126
 71   PTB(N,4)=SAVE(N)                                                  KPBP0127
 72   B3=PTB(7,3)*(PTB(3,3)-PTB(J+30,3))/PTB(J+30,3)                    KPBP0128
      B4=PTB(7,4)*(PTB(3,4)-PTB(J+30,4))/PTB(J+30,4)                    KPBP0129
      D3=PTB(7,3)*(PTB(4,3)-PTB(J+40,3))/PTB(J+30,3)                    KPBP0130
      D4=PTB(7,4)*(PTB(4,4)-PTB(J+40,4))/PTB(J+30,4)                    KPBP0131
      Z3Z4=PTB(6,3)-PTB(6,4)                                            KPBP0132
      PTB(J+30,3)=PTB(J+30,4)+A(J)*Z3Z4*(B3+B4)/2.                      KPBP0133
      PTB(J+40,3)=PTB(J+40,4)+A(J)*Z3Z4*(D3+D4)/2.                      KPBP0134
      PTB(J+10,3)=PTB(J+10,4)+Z3Z4*(EFN(J,3)+EFN(J,4))/2.               KPBP0135
      IF(J-JMAX) 801,800,800                                            KPBP0136
 800  PTB(J+20,3)=(PTB(J+20,4)*(PTB(J+30,4)*(PTB(5,4)**2-PTB(5,2)**2)   KPBP0137
     1   -PTB(J+40,4)*(PTB(6,4)-PTB(6,2))*(PTB(5,2)+PTB(5,4)))+         KPBP0138
     2   PTB(J+20,2)*(PTB(J+30,2)*(PTB(5,4)**2-PTB(5,3)**2)-            KPBP0139
     3   PTB(J+40,2)*((PTB(6,4)-PTB(6,2))*(PTB(5,4)+PTB(5,2))-          KPBP0140
     4   (PTB(6,3)-PTB(6,2))*(PTB(5,2)+PTB(5,3)))))/                    KPBP0141
     5   (PTB(J+30,3)*(PTB(5,3)**2-PTB(5,2)**2)-PTB(J+40,3)*            KPBP0142
     6   (PTB(6,3)-PTB(6,2))*(PTB(5,2)+PTB(5,3)))                       KPBP0143
      GO TO 89                                                          KPBP0144
 801  PTB(J+20,3)=RHO4(J)*(U4(J)*(PTB(5,1)**2-PTB(5,2)**2)              KPBP0145
     1   -V4(J)*((Z4-PTB(6,2))*(R4+PTB(5,2))+(PTB(6,1)-Z4)*             KPBP0146
     2   (R4+PTB(5,1))))+PTB(J+20,1)*(PTB(J+30,1)*(PTB(5,3)**2          KPBP0147
     3   -R4**2)-PTB(J+40,1)*((PTB(6,1)-Z4)*(R4+PTB(5,1))-              KPBP0148
     4   (PTB(6,1)-PTB(6,3))*(PTB(5,1)+PTB(5,3))))+PTB(J+20,2)*         KPBP0149
     5   (PTB(J+30,2)*(R4**2-PTB(5,3)**2)-PTB(J+40,2)*(                 KPBP0150
     6   (Z4-PTB(6,2))*(R4+PTB(5,2))-(PTB(6,3)-PTB(6,2))*               KPBP0151
     7   (PTB(5,2)+PTB(5,3))))                                          KPBP0152
      PTB(J+20,3)=PTB(J+20,3)/(PTB(J+30,3)*(PTB(5,1)**2-PTB(5,2)**2)    KPBP0153
     1   -PTB(J+40,3)*((PTB(6,3)-PTB(6,2))*(PTB(5,2)+PTB(5,3))+         KPBP0154
     2   (PTB(6,1)-PTB(6,3))*(PTB(5,1)+PTB(5,3))))                      KPBP0155
 89   CONTINUE                                                          KPBP0156
C                                                                       KPBP0157
C     STEP 9 CALCULATION OF THE POINT 5 BUCKET                          KPBP0158
C                                                                       KPBP0159
      CALL TAFN(PTB(4,3)/PTB(3,3),PTB(4,5)/PTB(3,5),BARC)               KPBP0160
      PTB(5,5)=(PTB(5,3)+BARC*(PTB(6,2)-PTB(6,3)-PTB(5,2)               KPBP0161
     X         *(PTB(6,1)-PTB(6,2))/(PTB(5,1)-PTB(5,2))))               KPBP0162
     X         /(1.-BARC*(PTB(6,1)-PTB(6,2))/(PTB(5,1)-PTB(5,2)))       KPBP0163
      Q=(PTB(5,5)-PTB(5,2))/(PTB(5,1)-PTB(5,2))                         KPBP0164
      PTB(6,5)=PTB(6,2)+Q*(PTB(6,1)-PTB(6,2))                           KPBP0165
      PTB(1,5)=PTB(1,2)+Q*(PTB(1,1)-PTB(1,2))                           KPBP0166
      PTB(2,5)=PTB(2,2)+Q*(PTB(2,1)-PTB(2,2))                           KPBP0167
      PTB(3,5)=PTB(3,2)+Q*(PTB(3,1)-PTB(3,2))                           KPBP0168
      PTB(4,5)=PTB(4,2)+Q*(PTB(4,1)-PTB(4,2))                           KPBP0169
      PTB(7,5)=(PTB(1,5)/PTB(2,5))**CAPN                                KPBP0170
      IF(L) 94,94,92                                                    KPBP0171
 92   DO 93 J=1,L                                                       KPBP0172
      PTB(J+10,5)=PTB(J+10,2)+Q*(PTB(J+10,1)-PTB(J+10,2))               KPBP0173
      PTB(J+20,5)=PTB(J+20,2)+Q*(PTB(J+20,1)-PTB(J+20,2))               KPBP0174
      PTB(J+30,5)=PTB(J+30,2)+Q*(PTB(J+30,1)-PTB(J+30,2))               KPBP0175
 93   PTB(J+40,5)=PTB(J+40,2)+Q*(PTB(J+40,1)-PTB(J+40,2))               KPBP0176
 94   FLG5=0.                                                           KPBP0177
      IF ((PTB(5,3)-PTB(5,4))/(PTB(6,3)-PTB(6,4))-                      KPBP0178
     X    (PTB(5,3)-PTB(5,5))/(PTB(6,3)-PTB(6,5)))   95,95,100          KPBP0179
 95   PTB(JMAX+10,5)=(PTB(JMAX+10,2)+PTB(JMAX+10,4))/2.                 KPBP0180
      PTB(JMAX+20,5)=(PTB(JMAX+20,2)+PTB(JMAX+20,4))/2.                 KPBP0181
      PTB(JMAX+30,5)=(PTB(JMAX+30,2)+PTB(JMAX+30,4))/2.                 KPBP0182
      PTB(JMAX+40,5)=(PTB(JMAX+40,2)+PTB(JMAX+40,4))/2.                 KPBP0183
      FLG5=1.                                                           KPBP0184
C                                                                       KPBP0185
C     STEP 10 CALCULATION OF THE GAS PROPERTIES P,RHO,U, AND V          KPBP0186
C                                                                       KPBP0187
 100  ETA1=PTB(2,1)*(PTB(4,1)*RT1(1)-PTB(3,1)*RT2(1))/(RT1(1)*VAL1(1))  KPBP0188
      C8=PTB(4,3)*RT1(3)                                                KPBP0189
      C9=PTB(3,3)*RT2(3)                                                KPBP0190
      C10=RT1(3)*VAL1(3)/PTB(2,3)                                       KPBP0191
      ETA3=(C8-C9)/C10                                                  KPBP0192
      GMU3=(C8+C9)/C10                                                  KPBP0193
      CALL SUMP2(1)                                                     KPBP0194
      CALL SUMP2(2)                                                     KPBP0195
      CALL SUMP2(3)                                                     KPBP0196
      CALL SUMP2(5)                                                     KPBP0197
      F1  =RT3(1)*(PTB(7,1)*TRL4(1)-GAMMA*PTB(1,1)*PTB(4,1)/PTB(5,1))   KPBP0198
      F3  =RT3(3)*(PTB(7,3)*TRL4(3)-GAMMA*PTB(1,3)*PTB(4,3)/PTB(5,3))   KPBP0199
      F2K =RT3(2)*(PTB(7,2)*TRM4(2)-GAMMA*PTB(1,2)*PTB(4,2)/PTB(5,2))   KPBP0200
      F3K =RT3(3)*(PTB(7,3)*TRM4(3)-GAMMA*PTB(1,3)*PTB(4,3)/PTB(5,3))   KPBP0201
      C11=PTB(7,3)*TRM1(3)                                              KPBP0202
      C12=PTB(7,3)*TRM2(3)                                              KPBP0203
      C13=PTB(7,3)*TRM3(3)                                              KPBP0204
      C19=PTB(7,3)*TRL1(3)                                              KPBP0205
      C20=PTB(7,3)*TRL2(3)                                              KPBP0206
      C21=PTB(7,3)*TRL3(3)                                              KPBP0207
      CPJK2=PTB(1,2)+(PTB(6,3)-PTB(6,2))*(F3K *GMU3+F2K *GMU2)          KPBP0208
     X      /2.+(RT3(2)+RT3(3))*(                                       KPBP0209
     X      PTB(4,2)*(PTB(2,2)*PTB(3,2)+PTB(2,3)*PTB(3,3))-             KPBP0210
     X      PTB(3,2)*(PTB(2,2)*PTB(4,2)+PTB(2,3)*PTB(4,3))              KPBP0211
     X      -(PTB(5,3)-PTB(5,2))*(C12+PTB(7,2)*TRM2(2))                 KPBP0212
     X      +(PTB(6,3)-PTB(6,2))*(C13+PTB(7,2)*TRM3(2))                 KPBP0213
     X      +(PTB(3,2)*(PTB(5,3)-PTB(5,2))-PTB(4,2)*(PTB(6,3)-PTB(6,2)) KPBP0214
     X      -PTB(3,3)*PTB(5,2)+PTB(4,3)*PTB(6,2))/2.*                   KPBP0215
     X      (C11+PTB(7,2)*TRM1(2)))/4.                                  KPBP0216
      CAPJ3=PTB(1,1)-(PTB(6,3)-PTB(6,1))*(F3  *ETA3+F1  *ETA1)          KPBP0217
     X      /2.-(RT3(1)+RT3(3))*(                                       KPBP0218
     X      PTB(4,1)*(PTB(2,1)*PTB(3,1)+PTB(2,3)*PTB(3,3))-             KPBP0219
     X      PTB(3,1)*(PTB(2,1)*PTB(4,1)+PTB(2,3)*PTB(4,3))              KPBP0220
     X      -(PTB(5,3)-PTB(5,1))*(C20+PTB(7,1)*TRL2(1))                 KPBP0221
     X      +(PTB(6,3)-PTB(6,1))*(C21+PTB(7,1)*TRL3(1))                 KPBP0222
     X      +(PTB(3,1)*(PTB(5,3)-PTB(5,1))-PTB(4,1)*(PTB(6,3)-PTB(6,1)) KPBP0223
     X      -PTB(3,3)*PTB(5,1)+PTB(4,3)*PTB(6,1))/2.*                   KPBP0224
     X      (C19+PTB(7,1)*TRL1(1)))/4.                                  KPBP0225
      CPLK2=RT3(2)*(PTB(2,2)*PTB(4,2)+PTB(5,3)*PTB(7,2)*TRM1(2)/2.)     KPBP0226
      CPLK3=RT3(3)*(PTB(2,3)*PTB(4,3)+PTB(5,3)*PTB(7,3)*TRM1(3)/2.)     KPBP0227
      CAPL(1)=RT3(1)*(PTB(2,1)*PTB(4,1)+PTB(5,3)*PTB(7,1)*TRL1(1)/2.)   KPBP0228
      CAPL(3)=RT3(3)*(PTB(2,3)*PTB(4,3)+PTB(5,3)*PTB(7,3)*TRL1(3)/2.)   KPBP0229
      QK2= -RT3(2)*(PTB(2,2)*PTB(3,2)+PTB(6,3)*PTB(7,2)*TRM1(2)/2.)     KPBP0230
      QK3= -RT3(3)*(PTB(2,3)*PTB(3,3)+PTB(6,3)*PTB(7,3)*TRM1(3)/2.)     KPBP0231
      Q(1)=  -RT3(1)*(PTB(2,1)*PTB(3,1)+PTB(6,3)*PTB(7,1)*TRL1(1)/2.)   KPBP0232
      Q(3)=  -RT3(3)*(PTB(2,3)*PTB(3,3)+PTB(6,3)*PTB(7,3)*TRL1(3)/2.)   KPBP0233
      S3=PTB(2,3)/(GAMMA*PTB(1,3))                                      KPBP0234
      S5=PTB(2,5)/(GAMMA*PTB(1,5))                                      KPBP0235
      IF(FLG5) 102,102,101                                              KPBP0236
 101  CAPJ1=PTB(1,5)+((PTB(5,3)-PTB(5,5))*(C13+PTB(7,5)*TRM3(5))        KPBP0237
     X               +(PTB(6,3)-PTB(6,5))*(C12+PTB(7,5)*TRM2(5))        KPBP0238
     X   -PTB(3,5)/2.*(PTB(6,3)-PTB(6,5))*(C11+PTB(7,5)*TRM1(5))        KPBP0239
     X   -PTB(4,5)/2.*(PTB(5,3)-PTB(5,5))*(C11+PTB(7,5)*TRM1(5))        KPBP0240
     X      +PTB(4,5)*(PTB(2,3)*PTB(4,3)+PTB(2,5)*PTB(4,5))             KPBP0241
     X      +PTB(3,5)*(PTB(2,3)*PTB(3,3)+PTB(2,5)*PTB(3,5)))/2.         KPBP0242
      G3=PTB(2,3)*PTB(3,3)+(PTB(6,3)-PTB(6,5))*C11/2.                   KPBP0243
      H3=PTB(2,3)*PTB(4,3)+(PTB(5,3)-PTB(5,5))*C11/2.                   KPBP0244
      G5=PTB(2,5)*PTB(3,5)+(PTB(6,3)-PTB(6,5))*PTB(7,5)*TRM1(5)/2.      KPBP0245
      H5=PTB(2,5)*PTB(4,5)+(PTB(5,3)-PTB(5,5))*PTB(7,5)*TRM1(5)/2.      KPBP0246
      T3=S3/PTB(3,3)*PTB(7,3)*TRM4(3)                                   KPBP0247
      T5=S5/PTB(3,5)*PTB(7,5)*TRM4(5)                                   KPBP0248
      GO TO 103                                                         KPBP0249
 102  CAPJ1=PTB(1,5)+((PTB(5,3)-PTB(5,5))*(C21+PTB(7,5)*TRL3(5))        KPBP0250
     X               +(PTB(6,3)-PTB(6,5))*(C20+PTB(7,5)*TRL2(5))        KPBP0251
     X   -PTB(3,5)/2.*(PTB(6,3)-PTB(6,5))*(C19+PTB(7,5)*TRL1(5))        KPBP0252
     X   -PTB(4,5)/2.*(PTB(5,3)-PTB(5,5))*(C19+PTB(7,5)*TRL1(5))        KPBP0253
     X      +PTB(4,5)*(PTB(2,3)*PTB(4,3)+PTB(2,5)*PTB(4,5))             KPBP0254
     X      +PTB(3,5)*(PTB(2,3)*PTB(3,3)+PTB(2,5)*PTB(3,5)))/2.         KPBP0255
      G3=PTB(2,3)*PTB(3,3)+(PTB(6,3)-PTB(6,5))*C19/2.                   KPBP0256
      H3=PTB(2,3)*PTB(4,3)+(PTB(5,3)-PTB(5,5))*C19/2.                   KPBP0257
      G5=PTB(2,5)*PTB(3,5)+(PTB(6,3)-PTB(6,5))*PTB(7,5)*TRL1(5)/2.      KPBP0258
      H5=PTB(2,5)*PTB(4,5)+(PTB(5,3)-PTB(5,5))*PTB(7,5)*TRL1(5)/2.      KPBP0259
      T3=S3/PTB(3,3)*PTB(7,3)*TRL4(3)                                   KPBP0260
      T5=S5/PTB(3,5)*PTB(7,5)*TRL4(5)                                   KPBP0261
 103  C14=CAPJ3-CPJK2                                                   KPBP0262
      C15=QK3+QK2+H5+H3                                                 KPBP0263
      C16=Q(3)+QK3+Q(1)+QK2                                             KPBP0264
      C17=CAPL(3)+CPLK3+CAPL(1)+CPLK2                                   KPBP0265
      C18=CPLK3+CPLK2+G5+G3                                             KPBP0266
      PTB(3,3)=2.*(C14*C15-(CAPJ1-CPJK2)*C16)/(C17*C15-C18*C16)         KPBP0267
      PTB(4,3)=(2.*C14-PTB(3,3)*C17)/C16                                KPBP0268
      PTB(1,3)=CPJK2+(PTB(3,3)*(CPLK3+CPLK2)+PTB(4,3)*                  KPBP0269
     X         (QK3+QK2))/2.                                            KPBP0270
      PTB(2,3)=PTB(2,5)+((S3+S5)*(PTB(1,3)-PTB(1,5))-(T3+T5)            KPBP0271
     X         *(PTB(6,3)-PTB(6,5)))/2.                                 KPBP0272
      CALL CHECK(DIFF)                                                  KPBP0273
      IF(DIFF) 80,81,80                                                 KPBP0274
 80   CONTINUE                                                          KPBP0275
      I=IMAX                                                            KPBP0276
 81   ICNT=I                                                            KPBP0277
      PTB(8,3)=1.                                                       KPBP0278
      RETURN                                                            KPBP0279
      END                                                               KPBP0280
$IBFTC PRINT                                                            PRIN0000
      SUBROUTINE PRINT (ID)                                             PRIN0001
C                                                                       PRIN0002
      COMMON/NAMEP/P(50,100)                                            PRIN0003
      DIMENSION   A(10)    ,BS(10)   ,WPWG(10)           ,PP(50,100),   PRIN0004
     X PS(50,50) ,PTB(50,8),PW(150)  ,R(10)    ,TRL1(5)  ,TRL2(5)   ,   PRIN0005
     X TRL3(5)   ,TRL4(5)  ,TRM1(5)  ,TRM2(5)  ,TRM3(5)  ,TRM4(5)   ,   PRIN0006
     X FA(100)   ,AAA(10)  ,WDP(10)  ,AA(8)                             PRIN0007
     X,REK(10)   ,VPK(10)  ,THA(10)  ,TPK(10)  ,RPRG(10) ,RPRP0(10)     PRIN0008
     X,XXX(159)                                                         PRIN0009
     X,SP(10)                                                           PRIN0010
      COMMON CAPN  ,CPL   ,CPS   ,CTOP  ,DL    ,DR    ,INSERT,EKW   ,   PRIN0011
     XEW    ,WPWG  ,GAMMA ,GMG0  ,HPL   ,HPS   ,IMAX  ,N1    ,N2    ,   PRIN0012
     XPS    ,PG0   ,PW    ,R     ,RCAP  ,RHOP0 ,SMP   ,TG0   ,TPM   ,   PRIN0013
     XA     ,BS    ,CC    ,CONST ,EOW1  ,EOW2  ,HIT   ,ICNT  ,IPP   ,   PRIN0014
     XJMAX  ,MW           ,PP    ,PTB   ,TRL1  ,TRL2  ,TRL3  ,TRL4  ,   PRIN0015
     XTRM1  ,TRM2  ,TRM3  ,TRM4  ,LRC   ,DELTA ,FA    ,VGR   ,CPG   ,   PRIN0016
     XNILP  ,THI   ,AAA   ,WDG   ,IFG   ,RHOG0 ,THRUST,TI    ,WDP   ,   PRIN0017
     XAA    ,RT    ,VK    ,T     ,XXX   ,NPLOT                          PRIN0018
C                                                                       PRIN0019
      I=3                                                               PRIN0020
      II=0                                                              PRIN0021
      IE=ID                                                             PRIN0022
      ICN=ICNT                                                          PRIN0023
      CF=0.                                                             PRIN0024
      XISP=0.                                                           PRIN0025
      GO TO (1,2,3,3,5,6,6),ID                                          PRIN0026
 1    ICN=0                                                             PRIN0027
      IF(LRC) 12,11,12                                                  PRIN0028
 11   I=1                                                               PRIN0029
      II=14                                                             PRIN0030
      LRC1=1                                                            PRIN0031
      WDPS=0.                                                           PRIN0032
      DO 111 M=1,10                                                     PRIN0033
      IF(WDP(M)) 111,5,111                                              PRIN0034
 111  WDPS=WDPS+WDP(M)                                                  PRIN0035
      GO TO 5                                                           PRIN0036
 12   I=2                                                               PRIN0037
      II=12                                                             PRIN0038
      IF(PTB(5,2)) 103,3,103                                            PRIN0039
 103  II=13                                                             PRIN0040
      IF(PTB(8,2)) 3,30,3                                               PRIN0041
 2    N3=N3+1                                                           PRIN0042
      IF( MOD (N3,N2)) 999,203,999                                      PRIN0043
 203  IF( MOD (LRC,N1)) 999,30,999                                      PRIN0044
 5    CF=THRUST*10.241302/RT/RT/PG0                                     PRIN0045
      XISP=THRUST/(WDG+WDPS)                                            PRIN0046
      N3=0                                                              PRIN0047
      GO TO 3                                                           PRIN0048
 6    IF(N1+N2-2) 601,601,999                                           PRIN0049
 601  ICN=0                                                             PRIN0050
      I=ID-5                                                            PRIN0051
 3    IE=ID+II                                                          PRIN0052
 30   TG=PTB(1,I)/PTB(2,I)/RCAP                                         PRIN0053
      VG=SQRT (PTB(3,I)**2+PTB(4,I)**2)                                 PRIN0054
      ACHM=VG/SQRT (GAMMA*RCAP*TG)                                      PRIN0055
      THETA=57.29578*ATAN2(PTB(4,I),PTB(3,I))                           PRIN0056
      TGTG0=TG/TG0                                                      PRIN0057
      PGPG0=PTB(1,I)/PG0                                                PRIN0058
      RGRG0=PTB(2,I)/RHOG0                                              PRIN0059
      SRKRG=0.                                                          PRIN0060
      IF(JMAX) 39,39,31                                                 PRIN0061
 31   DO 32 J=1,JMAX                                                    PRIN0062
 32   SRKRG=SRKRG+PTB(J+20,I)                                           PRIN0063
      SRKRG=SRKRG/PTB(2,I)                                              PRIN0064
      DO 35 N=1,JMAX                                                    PRIN0065
      REK(N)=2.*PTB(2,I)*R(N)*SQRT ((PTB(3,I)-PTB(N+30,I))              PRIN0066
     X**2 +(PTB(4,I)-PTB(N+40,I))**2 )/PTB(7,I)/CC                      PRIN0067
      VPK(N)=SQRT (PTB(N+30,I)**2 +PTB(N+40,I)**2)                      PRIN0068
      THA(N)=57.29578*ATAN2(PTB(N+40,I),PTB(N+30,I))                    PRIN0069
      IF (PTB(N+10,I)-HPL) 300,300,303                                  PRIN0070
 300  IF (PTB(N+10,I)-HPS) 301,302,302                                  PRIN0071
 301  TPK(N)=PTB(N+10,I)/CPS                                            PRIN0072
      GO TO 304                                                         PRIN0073
 302  TPK(N)=TPM                                                        PRIN0074
      GO TO 304                                                         PRIN0075
 303  TPK(N)=TPM+(PTB(N+10,I)-HPL)/CPL                                  PRIN0076
 304  RPRG(N)=PTB(N+20,I)/PTB(2,I)                                      PRIN0077
 35   RPRP0(N)=PTB(N+20,I)/RHOP0                                        PRIN0078
 39   IF(LRC-LRC1) 36,37,36                                             PRIN0079
 36   WRITE (6,98)                                                      PRIN0080
      WRITE (6,99)                                                      PRIN0081
      LRC1=LRC                                                          PRIN0082
 37   WRITE (6,100)LRC,IE,PTB(5,I),PTB(6,I),ACHM,                       PRIN0083
     XTG,VG,THETA,TGTG0,PGPG0,RGRG0,SRKRG,CF,XISP,ICN                   PRIN0084
      IF(JMAX) 999,999,38                                               PRIN0085
 38   WRITE (6,101)(L,REK(L),VPK(L),THA(L),                             PRIN0086
     XTPK(L),RPRG(L),RPRP0(L),R(L),L=1,JMAX)                            PRIN0087
 98   FORMAT (7H0LRC ID,4X,1HR,8X,1HZ,7X,4HMACH,5X,2HTG,5X,2HVG,        PRIN0088
     X 4X,7HTHETA-G,2X,6HTG/TG0,2X,6HPG/PG0,2X,6HDG/DG0,2X,             PRIN0089
     X 6HSDK/DG,3X,2HCF,6X,3HISP,3X,2HIT)                               PRIN0090
 99   FORMAT (8X,1HK,6X,3HREK,11X,3HVPK,10X,7HTHETA-K,8X,               PRIN0091
     X 3HTPK,8X,6HDPK/DG,8X,7HDPK/DP0,10X,3HRPK)                        PRIN0092
 100  FORMAT (1H0,I3,I3, F9.5, F9.5,F8.4,F8.1 ,F8.1 ,F8.3,F8.5,F8.5,F8.5PRIN0093
     X ,F8.5,F8.4, F8.2,I3)                                             PRIN0094
 101  FORMAT (I9,7E14.6)                                                PRIN0095
C                                                                       PRIN0096
 999  RETURN                                                            PRIN0097
       END                                                              PRIN0098
$IBFTC PTINT                                                            PTIN0000
      SUBROUTINE PTINT(IF)                                              PTIN0001
C     THE INTERIOR POINT CALCULATION                                    PTIN0002
C                                                                       PTIN0003
      COMMON/NAMEP/P(50,100)                                            PTIN0004
      DIMENSION   A(10)    ,BS(10)   ,FJ(10)             ,PP(50,100),   PTIN0005
     X PS(50,50) ,PTB(50,8),PW(150)  ,R(10)    ,TRL1(5)  ,TRL2(5)   ,   PTIN0006
     X TRL3(5)   ,TRL4(5)  ,TRM1(5)  ,TRM2(5)  ,TRM3(5)  ,TRM4(5)   ,   PTIN0007
     XFA(100)                                                           PTIN0008
     X,RT1(3)    ,RT2(3)   ,RT3(3)   ,VAL1(3)  ,CAPL(3)                 PTIN0009
     X,Q(3)                                                             PTIN0010
     X,RHO4(10),U4(10),V4(10)                                           PTIN0011
      COMMON CAPN  ,CPL   ,CPS   ,CTOP  ,DL    ,DR    ,EK    ,EKW   ,   PTIN0012
     XEW    ,FJ    ,GAMMA ,GMG0  ,HPL   ,HPS   ,IMAX  ,N1    ,N2    ,   PTIN0013
     XPS    ,PG0   ,PW    ,R     ,RCAP  ,RHOP0 ,SMP   ,TG0   ,TPM   ,   PTIN0014
     XA     ,BS    ,CC    ,CONST ,EOW1  ,EOW2  ,HIT   ,ICNT  ,IPP   ,   PTIN0015
     XJMAX  ,MW           ,PP    ,PTB   ,TRL1  ,TRL2  ,TRL3  ,TRL4  ,   PTIN0016
     XTRM1  ,TRM2  ,TRM3  ,TRM4  ,LRC   ,DELTA ,FA    ,VGR              PTIN0017
C                                                                       PTIN0018
C     STEP 1 AVERAGEING FOR INITIAL CONDITIONS                          PTIN0019
C                                                                       PTIN0020
      DO 1 I=1,4                                                        PTIN0021
 1    PTB(I,3)=(PTB(I,1)+PTB(I,2))/2.                                   PTIN0022
      PTB(3,5)=PTB(3,3)                                                 PTIN0023
      PTB(4,5)=PTB(4,3)                                                 PTIN0024
      DO 11 N=1,2                                                       PTIN0025
      RT1(N)=SQRT (GAMMA*PTB(1,N))                                      PTIN0026
      RT2(N)=SQRT (PTB(2,N)*(PTB(3,N)**2 +PTB(4,N)**2 )-GAMMA*PTB(1,N)) PTIN0027
      RT3(N)=RT1(N)/RT2(N)                                              PTIN0028
 11   VAL1(N)=PTB(3,N)**2 *PTB(2,N)-GAMMA*PTB(1,N)                      PTIN0029
      GK2=(PTB(2,2)*PTB(3,2)*PTB(4,2)+RT1(2)*RT2(2))/VAL1(2)            PTIN0030
      GL1=(PTB(2,1)*PTB(3,1)*PTB(4,1)-RT1(1)*RT2(1))/VAL1(1)            PTIN0031
      ETA1=PTB(2,1)*(PTB(4,1)*RT1(1)-PTB(3,1)*RT2(1))/(RT1(1)*VAL1(1))  PTIN0032
      GMU2=PTB(2,2)*(PTB(4,2)*RT1(2)+PTB(3,2)*RT2(2))/(RT1(2)*VAL1(2))  PTIN0033
      IF(JMAX) 2,2,12                                                   PTIN0034
 12   DO 13 J=1,JMAX                                                    PTIN0035
      RHO4(J)=PTB(J+20,4)                                               PTIN0036
      U4(J)=PTB(J+30,4)                                                 PTIN0037
      V4(J)=PTB(J+40,4)                                                 PTIN0038
      PTB(J+10,3)=(PTB(J+10,1)+PTB(J+10,2))/2.                          PTIN0039
      PTB(J+20,3)=(PTB(J+20,1)+PTB(J+20,2))/2.                          PTIN0040
      PTB(J+30,3)=(PTB(J+30,1)+PTB(J+30,2))/2.                          PTIN0041
      PTB(J+40,3)=(PTB(J+40,1)+PTB(J+40,2))/2.                          PTIN0042
      PTB(J+30,4)=PTB(J+30,3)                                           PTIN0043
 13   PTB(J+40,4)=PTB(J+40,3)                                           PTIN0044
      R4=PTB(5,4)                                                       PTIN0045
      Z4=PTB(6,4)                                                       PTIN0046
      CALL ACOMP(1,2)                                                   PTIN0047
      CALL SUMP1(1)                                                     PTIN0048
      CALL SUMP1(2)                                                     PTIN0049
      F1  =RT3(1)*(PTB(7,1)*TRM4(1)-GAMMA*PTB(1,1)*PTB(4,1)/PTB(5,1))   PTIN0050
      IF(PTB(5,2)) 10,101,10                                            PTIN0051
 101  F2  =RT3(2)*(PTB(7,2)*TRM4(2)-GAMMA*PTB(1,2)*VGR)                 PTIN0052
      GO TO 3                                                           PTIN0053
 10   F2  =RT3(2)*(PTB(7,2)*TRM4(2)-GAMMA*PTB(1,2)*PTB(4,2)/PTB(5,2))   PTIN0054
      GO TO 3                                                           PTIN0055
C                                                                       PTIN0056
 2    CAPL(1)=RT3(1)*PTB(2,1)*PTB(4,1)                                  PTIN0057
      CAPL(2)=RT3(2)*PTB(2,2)*PTB(4,2)                                  PTIN0058
      Q(2)=-RT3(2)*PTB(2,2)*PTB(3,2)                                    PTIN0059
      F1  =-RT3(1)*GAMMA*PTB(1,1)*PTB(4,1)/PTB(5,1)                     PTIN0060
      IF(PTB(5,2)) 20,201,20                                            PTIN0061
 20   F2  =-RT3(2)*GAMMA*PTB(1,2)*VGR                                   PTIN0062
      GO TO 3                                                           PTIN0063
 201  F2  =-RT3(2)*GAMMA*PTB(1,2)*PTB(4,2)/PTB(5,2)                     PTIN0064
C                                                                       PTIN0065
C     STEP 2 CALCULATION OF BARK, BARL, Z3,AND R3                       PTIN0066
C            BEGIN ITERATION LOOP                                       PTIN0067
C                                                                       PTIN0068
 3    DO 80 I=1,IMAX                                                    PTIN0069
      DO 21 N=1,4                                                       PTIN0070
 21   PTB(N,8)=PTB(N,3)                                                 PTIN0071
      RT1(3)=SQRT (GAMMA*PTB(1,3))                                      PTIN0072
      RT2(3)=SQRT (PTB(2,3)*(PTB(3,3)**2 +PTB(4,3)**2 )-GAMMA*PTB(1,3)) PTIN0073
      RT3(3)=RT1(3)/RT2(3)                                              PTIN0074
      VAL1(3)=PTB(3,3)**2 *PTB(2,3)-GAMMA*PTB(1,3)                      PTIN0075
      IF(JMAX) 23,23,211                                                PTIN0076
 211  DO 22 J=1,JMAX                                                    PTIN0077
      PTB(J+10,8)=PTB(J+10,3)                                           PTIN0078
      PTB(J+20,8)=PTB(J+20,3)                                           PTIN0079
      PTB(J+30,8)=PTB(J+30,3)                                           PTIN0080
 22   PTB(J+40,8)=PTB(J+40,3)                                           PTIN0081
      PTB(7,3)=(PTB(1,3)/PTB(2,3))**CAPN                                PTIN0082
 23   GK3=(PTB(2,3)*PTB(3,3)*PTB(4,3)+RT1(3)*RT2(3))/VAL1(3)            PTIN0083
      CALL TAFN(GK2,GK3,BARK)                                           PTIN0084
      GL3=(PTB(2,3)*PTB(3,3)*PTB(4,3)-RT1(3)*RT2(3))/VAL1(3)            PTIN0085
      CALL TAFN(GL1,GL3,BARL)                                           PTIN0086
      PTB(6,3)=(PTB(5,1)-PTB(5,2)+BARK*PTB(6,2)-BARL*PTB(6,1))          PTIN0087
     X         /(BARK-BARL)                                             PTIN0088
      PTB(5,3)=PTB(5,2)+BARK*(PTB(6,3)-PTB(6,2))                        PTIN0089
C                                                                       PTIN0090
C     STEP 3 CALCULATION OF THE POINT 5 BUCKET                          PTIN0091
C                                                                       PTIN0092
      CALL TAFN(PTB(4,3)/PTB(3,3),PTB(4,5)/PTB(3,5),BARC)               PTIN0093
      PTB(5,5)=(PTB(5,3)+BARC*(PTB(6,2)-PTB(6,3)-PTB(5,2)               PTIN0094
     X         *(PTB(6,1)-PTB(6,2))/(PTB(5,1)-PTB(5,2))))               PTIN0095
     X         /(1.-BARC*(PTB(6,1)-PTB(6,2))/(PTB(5,1)-PTB(5,2)))       PTIN0096
      Q=(PTB(5,5)-PTB(5,2))/(PTB(5,1)-PTB(5,2))                         PTIN0097
      PTB(6,5)=PTB(6,2)+Q*(PTB(6,1)-PTB(6,2))                           PTIN0098
      PTB(1,5)=PTB(1,2)+Q*(PTB(1,1)-PTB(1,2))                           PTIN0099
      PTB(2,5)=PTB(2,2)+Q*(PTB(2,1)-PTB(2,2))                           PTIN0100
      PTB(3,5)=PTB(3,2)+Q*(PTB(3,1)-PTB(3,2))                           PTIN0101
      PTB(4,5)=PTB(4,2)+Q*(PTB(4,1)-PTB(4,2))                           PTIN0102
      C8=PTB(4,3)*RT1(3)                                                PTIN0103
      C9=PTB(3,3)*RT2(3)                                                PTIN0104
      C10=RT1(3)*VAL1(3)/PTB(2,3)                                       PTIN0105
      ETA3=(C8-C9)/C10                                                  PTIN0106
      GMU3=(C8+C9)/C10                                                  PTIN0107
      IF(JMAX) 75,75,71                                                 PTIN0108
 71   PTB(7,5)=(PTB(1,5)/PTB(2,5))**CAPN                                PTIN0109
      DO 69 J=1,JMAX                                                    PTIN0110
      PTB(J+10,5)=PTB(J+10,2)+Q*(PTB(J+10,1)-PTB(J+10,2))               PTIN0111
      PTB(J+20,5)=PTB(J+20,2)+Q*(PTB(J+20,1)-PTB(J+20,2))               PTIN0112
      PTB(J+30,5)=PTB(J+30,2)+Q*(PTB(J+30,1)-PTB(J+30,2))               PTIN0113
      PTB(J+40,5)=PTB(J+40,2)+Q*(PTB(J+40,1)-PTB(J+40,2))               PTIN0114
C                                                                       PTIN0115
C     STEP 4 CALCULATION OF THE POINT 4J BUCKET                         PTIN0116
C                                                                       PTIN0117
      CALL TAFN(PTB(J+40,3)/PTB(J+30,3),PTB(J+40,4)/PTB(J+30,4),BARC)   PTIN0118
      PTB(5,4)=(PTB(5,3)+BARC*(PTB(6,2)-PTB(6,3)-PTB(5,2)               PTIN0119
     X         *(PTB(6,1)-PTB(6,2))/(PTB(5,1)-PTB(5,2))))               PTIN0120
     X         /(1.-BARC*(PTB(6,1)-PTB(6,2))/(PTB(5,1)-PTB(5,2)))       PTIN0121
      Q=(PTB(5,4)-PTB(5,2))/(PTB(5,1)-PTB(5,2))                         PTIN0122
      PTB(6,4)=PTB(6,2)+Q*(PTB(6,1)-PTB(6,2))                           PTIN0123
      PTB(1,4)=PTB(1,2)+Q*(PTB(1,1)-PTB(1,2))                           PTIN0124
      PTB(2,4)=PTB(2,2)+Q*(PTB(2,1)-PTB(2,2))                           PTIN0125
      PTB(3,4)=PTB(3,2)+Q*(PTB(3,1)-PTB(3,2))                           PTIN0126
      PTB(4,4)=PTB(4,2)+Q*(PTB(4,1)-PTB(4,2))                           PTIN0127
      PTB(7,4)=(PTB(1,4)/PTB(2,4))**CAPN                                PTIN0128
      PTB(J+10,4)=PTB(J+10,2)+Q*(PTB(J+10,1)-PTB(J+10,2))               PTIN0129
      PTB(J+20,4)=PTB(J+20,2)+Q*(PTB(J+20,1)-PTB(J+20,2))               PTIN0130
      PTB(J+30,4)=PTB(J+30,2)+Q*(PTB(J+30,1)-PTB(J+30,2))               PTIN0131
      PTB(J+40,4)=PTB(J+40,2)+Q*(PTB(J+40,1)-PTB(J+40,2))               PTIN0132
C                                                                       PTIN0133
C     STEP 5 CALCULATION OF PARTICLE H,U,AND V AT PT. 3                 PTIN0134
C                                                                       PTIN0135
      B3=PTB(7,3)*(PTB(3,3)-PTB(J+30,3))/PTB(J+30,3)                    PTIN0136
      B4=PTB(7,4)*(PTB(3,4)-PTB(J+30,4))/PTB(J+30,4)                    PTIN0137
      D3=PTB(7,3)*(PTB(4,3)-PTB(J+40,3))/PTB(J+30,3)                    PTIN0138
      D4=PTB(7,4)*(PTB(4,4)-PTB(J+40,4))/PTB(J+30,4)                    PTIN0139
      Z3Z4=PTB(6,3)-PTB(6,4)                                            PTIN0140
      PTB(J+30,3)=PTB(J+30,4)+A(J)*Z3Z4*(B3+B4)/2.                      PTIN0141
      PTB(J+40,3)=PTB(J+40,4)+A(J)*Z3Z4*(D3+D4)/2.                      PTIN0142
      PTB(J+10,3)=PTB(J+10,4)+Z3Z4*(EFN(J,3)+EFN(J,4))/2.               PTIN0143
C                                                                       PTIN0144
C     STEP 6 CALCULATION OF PARTICLE RHO                                PTIN0145
C                                                                       PTIN0146
      IF(IF) 60,60,65                                                   PTIN0147
 60   PTB(J+20,3)=(PTB(J+20,2)*(PTB(J+30,2)*(PTB(5,1)**2-               PTIN0148
     X   PTB(5,3)**2)-PTB(J+40,2)*((PTB(6,1)-PTB(6,2))*(PTB(5,1)        PTIN0149
     X   +PTB(5,2))-(PTB(5,3)**2-PTB(5,2)**2)/GK2))+                    PTIN0150
     X   PTB(J+20,1)*(PTB(J+30,1)*(PTB(5,3)**2-PTB(5,2)**2)             PTIN0151
     X   -PTB(J+40,1)*((PTB(6,1)-PTB(6,2))*(PTB(5,1)+PTB(5,2))-         PTIN0152
     X   (PTB(5,1)+PTB(5,3))*(PTB(6,1)-PTB(6,3)))))/                    PTIN0153
     X   (PTB(J+30,3)*(PTB(5,1)**2-PTB(5,2)**2)-PTB(J+40,3)*            PTIN0154
     X   ((PTB(5,3)**2-PTB(5,2)**2)/GK3+                                PTIN0155
     X   (PTB(5,1)+PTB(5,3))*(PTB(6,1)-PTB(6,3))))                      PTIN0156
      GO TO 69                                                          PTIN0157
 65   PTB(J+20,3)=RHO4(J)*(U4(J)*(PTB(5,1)**2-PTB(5,2)**2)              PTIN0158
     1   -V4(J)*((Z4-PTB(6,2))*(R4+PTB(5,2))+(PTB(6,1)-Z4)*             PTIN0159
     2   (R4+PTB(5,1))))+PTB(J+20,1)*(PTB(J+30,1)*(PTB(5,3)**2          PTIN0160
     3   -R4**2)-PTB(J+40,1)*((PTB(6,1)-Z4)*(R4+PTB(5,1))-              PTIN0161
     4   (PTB(6,1)-PTB(6,3))*(PTB(5,1)+PTB(5,3))))+PTB(J+20,2)*         PTIN0162
     5   (PTB(J+30,2)*(R4**2-PTB(5,3)**2)-PTB(J+40,2)*(                 PTIN0163
     6   (Z4-PTB(6,2))*(R4+PTB(5,2))-(PTB(6,3)-PTB(6,2))*               PTIN0164
     7   (PTB(5,2)+PTB(5,3))))                                          PTIN0165
      PTB(J+20,3)=PTB(J+20,3)/(PTB(J+30,3)*(PTB(5,1)**2-PTB(5,2)**2)    PTIN0166
     1   -PTB(J+40,3)*((PTB(6,3)-PTB(6,2))*(PTB(5,2)+PTB(5,3))+         PTIN0167
     2   (PTB(6,1)-PTB(6,3))*(PTB(5,1)+PTB(5,3))))                      PTIN0168
 69   CONTINUE                                                          PTIN0169
C                                                                       PTIN0170
C     STEP 7 CALCULATION OF GAS PROPERTIES P, RHO, U, AND V             PTIN0171
C                                                                       PTIN0172
      CALL SUMP1(3)                                                     PTIN0173
      CALL SUMP1(5)                                                     PTIN0174
      F3  =RT3(3)*(PTB(7,3)*TRM4(3)-GAMMA*PTB(1,3)*PTB(4,3)/PTB(5,3))   PTIN0175
      C11=PTB(7,3)*TRM1(3)                                              PTIN0176
      C12=PTB(7,3)*TRM2(3)                                              PTIN0177
      C13=PTB(7,3)*TRM3(3)                                              PTIN0178
      CAPJ2=PTB(1,2)+(PTB(6,3)-PTB(6,2))*(F3  *GMU3+F2  *GMU2)          PTIN0179
     X      /2.+(RT3(2)+RT3(3))*(                                       PTIN0180
     X      PTB(4,2)*(PTB(2,2)*PTB(3,2)+PTB(2,3)*PTB(3,3))-             PTIN0181
     X      PTB(3,2)*(PTB(2,2)*PTB(4,2)+PTB(2,3)*PTB(4,3))              PTIN0182
     X      -(PTB(5,3)-PTB(5,2))*(C12+PTB(7,2)*TRM2(2))                 PTIN0183
     X      +(PTB(6,3)-PTB(6,2))*(C13+PTB(7,2)*TRM3(2))                 PTIN0184
     X      +(PTB(3,2)*(PTB(5,3)-PTB(5,2))-PTB(4,2)*(PTB(6,3)-PTB(6,2)) PTIN0185
     X      -PTB(3,3)*PTB(5,2)+PTB(4,3)*PTB(6,2))/2.*                   PTIN0186
     X      (C11+PTB(7,2)*TRM1(2)))/4.                                  PTIN0187
      CAPJ3=PTB(1,1)-(PTB(6,3)-PTB(6,1))*(F3  *ETA3+F1  *ETA1)          PTIN0188
     X      /2.-(RT3(1)+RT3(3))*(                                       PTIN0189
     X      PTB(4,1)*(PTB(2,1)*PTB(3,1)+PTB(2,3)*PTB(3,3))-             PTIN0190
     X      PTB(3,1)*(PTB(2,1)*PTB(4,1)+PTB(2,3)*PTB(4,3))              PTIN0191
     X      -(PTB(5,3)-PTB(5,1))*(C12+PTB(7,1)*TRM2(1))                 PTIN0192
     X      +(PTB(6,3)-PTB(6,1))*(C13+PTB(7,1)*TRM3(1))                 PTIN0193
     X      +(PTB(3,1)*(PTB(5,3)-PTB(5,1))-PTB(4,1)*(PTB(6,3)-PTB(6,1)) PTIN0194
     X      -PTB(3,3)*PTB(5,1)+PTB(4,3)*PTB(6,1))/2.*                   PTIN0195
     X      (C11+PTB(7,1)*TRM1(1)))/4.                                  PTIN0196
      CAPJ1=PTB(1,5)+((PTB(5,3)-PTB(5,5))*(C13+PTB(7,5)*TRM3(5))        PTIN0197
     X               +(PTB(6,3)-PTB(6,5))*(C12+PTB(7,5)*TRM2(5))        PTIN0198
     X   -PTB(3,5)/2.*(PTB(6,3)-PTB(6,5))*(C11+PTB(7,5)*TRM1(5))        PTIN0199
     X   -PTB(4,5)/2.*(PTB(5,3)-PTB(5,5))*(C11+PTB(7,5)*TRM1(5))        PTIN0200
     X      +PTB(4,5)*(PTB(2,3)*PTB(4,3)+PTB(2,5)*PTB(4,5))             PTIN0201
     X      +PTB(3,5)*(PTB(2,3)*PTB(3,3)+PTB(2,5)*PTB(3,5)))/2.         PTIN0202
      G3=PTB(2,3)*PTB(3,3)+(PTB(6,3)-PTB(6,5))*C11/2.                   PTIN0203
      G5=PTB(2,5)*PTB(3,5)+(PTB(6,3)-PTB(6,5))*PTB(7,5)*TRM1(5)/2.      PTIN0204
      H3=PTB(2,3)*PTB(4,3)+(PTB(5,3)-PTB(5,5))*C11/2.                   PTIN0205
      H5=PTB(2,5)*PTB(4,5)+(PTB(5,3)-PTB(5,5))*PTB(7,5)*TRM1(5)/2.      PTIN0206
      DO 73 N=1,3                                                       PTIN0207
      CAPL(N)=RT3(N)*(PTB(2,N)*PTB(4,N)+PTB(5,3)*PTB(7,N)*TRM1(N)/2.)   PTIN0208
 73   Q(N)=  -RT3(N)*(PTB(2,N)*PTB(3,N)+PTB(6,3)*PTB(7,N)*TRM1(N)/2.)   PTIN0209
      S3=PTB(2,3)/(GAMMA*PTB(1,3))                                      PTIN0210
      S5=PTB(2,5)/(GAMMA*PTB(1,5))                                      PTIN0211
      T3=S3/PTB(3,3)*PTB(7,3)*TRM4(3)                                   PTIN0212
      T5=S5/PTB(3,5)*PTB(7,5)*TRM4(5)                                   PTIN0213
      GO TO 79                                                          PTIN0214
 75   F3  =-RT3(3)*GAMMA*PTB(1,3)*PTB(4,3)/PTB(5,3)                     PTIN0215
      CAPJ2=PTB(1,2)+(PTB(6,3)-PTB(6,2))*(F3  *GMU3+F2  *GMU2)          PTIN0216
     X      /2.+(RT3(2)+RT3(3))*(                                       PTIN0217
     X      PTB(4,2)*(PTB(2,2)*PTB(3,2)+PTB(2,3)*PTB(3,3))-             PTIN0218
     X      PTB(3,2)*(PTB(2,2)*PTB(4,2)+PTB(2,3)*PTB(4,3)))/4.          PTIN0219
      CAPJ3=PTB(1,1)-(PTB(6,3)-PTB(6,1))*(F3  *ETA3+F1  *ETA1)          PTIN0220
     X      /2.-(RT3(1)+RT3(3))*(                                       PTIN0221
     X      PTB(4,1)*(PTB(2,1)*PTB(3,1)+PTB(2,3)*PTB(3,3))-             PTIN0222
     X      PTB(3,1)*(PTB(2,1)*PTB(4,1)+PTB(2,3)*PTB(4,3)))/4.          PTIN0223
      CAPJ1=PTB(1,5)+(PTB(4,5)*(PTB(2,3)*PTB(4,3)+PTB(2,5)*             PTIN0224
     X    PTB(4,5))+PTB(3,5)*(PTB(2,3)*PTB(3,3)+PTB(2,5)*PTB(3,5)))/2.  PTIN0225
      G3=PTB(2,3)*PTB(3,3)                                              PTIN0226
      G5=PTB(2,5)*PTB(3,5)                                              PTIN0227
      H3=PTB(2,3)*PTB(4,3)                                              PTIN0228
      H5=PTB(2,5)*PTB(4,5)                                              PTIN0229
      CAPL(3)=RT3(3)*PTB(2,3)*PTB(4,3)                                  PTIN0230
      Q(1)=-RT3(1)*PTB(2,1)*PTB(3,1)                                    PTIN0231
      Q(3)=-RT3(3)*PTB(2,3)*PTB(3,3)                                    PTIN0232
      S3=PTB(2,3)/(GAMMA*PTB(1,3))                                      PTIN0233
      S5=PTB(2,5)/(GAMMA*PTB(1,5))                                      PTIN0234
      T3=0.                                                             PTIN0235
      T5=0.                                                             PTIN0236
 79   C14=CAPJ3-CAPJ2                                                   PTIN0237
      C15=Q(3)+Q(2)+H5+H3                                               PTIN0238
      C16=2.*Q(3)+Q(1)+Q(2)                                             PTIN0239
      C17=2.*CAPL(3)+CAPL(1)+CAPL(2)                                    PTIN0240
      PTB(3,3)=2.*(C14*C15-(CAPJ1-CAPJ2)*C16)/                          PTIN0241
     X         (C17*C15-(CAPL(3)+CAPL(2)+G5+G3)*C16)                    PTIN0242
      PTB(4,3)=(2.*C14-PTB(3,3)*C17)/C16                                PTIN0243
      IF(PTB(4,3)) 791,792,792                                          PTIN0244
 791  PTB(4,3)=PTB(4,2)+(PTB(5,3)-PTB(5,2))/(PTB(5,1)-PTB(5,2))         PTIN0245
     X         *(PTB(4,1)-PTB(4,2))                                     PTIN0246
 792  PTB(1,3)=CAPJ2+(PTB(3,3)*(CAPL(3)+CAPL(2))+                       PTIN0247
     X         PTB(4,3)*(Q(3)+Q(2)))/2.                                 PTIN0248
      PTB(2,3)=PTB(2,5)+((S3+S5)*(PTB(1,3)-PTB(1,5))-(T3+T5)            PTIN0249
     X         *(PTB(6,3)-PTB(6,5)))/2.                                 PTIN0250
      CALL CHECK(DIFF)                                                  PTIN0251
      IF(DIFF) 80,81,80                                                 PTIN0252
 80   CONTINUE                                                          PTIN0253
      I=IMAX                                                            PTIN0254
 81   ICNT=I                                                            PTIN0255
      RETURN                                                            PTIN0256
       END                                                              PTIN0257
$IBFTC S1UMP                                                            S1UM0000
      SUBROUTINE SUMP1(N)                                               S1UM0001
C                                                                       S1UM0002
      COMMON/NAMEP/P(50,100)                                            S1UM0003
      DIMENSION   A(10)    ,BS(10)   ,WPWG(10)           ,PP(50,100),   S1UM0004
     X PS(50,50) ,PTB(50,8),PW(150)  ,R(10)    ,TRL1(5)  ,TRL2(5)   ,   S1UM0005
     X TRL3(5)   ,TRL4(5)  ,TRM1(5)  ,TRM2(5)  ,TRM3(5)  ,TRM4(5)   ,   S1UM0006
     X FA(100)                                                          S1UM0007
      COMMON CAPN  ,CPL   ,CPS   ,CTOP  ,DL    ,DR    ,EK    ,EKW   ,   S1UM0008
     XEW    ,WPWG  ,GAMMA ,GMG0  ,HPL   ,HPS   ,IMAX  ,N1    ,N2    ,   S1UM0009
     XPS    ,PG0   ,PW    ,R     ,RCAP  ,RHOP0 ,SMP   ,TG0   ,TPM   ,   S1UM0010
     XA     ,BS    ,CC    ,CONST ,EOW1  ,EOW2  ,HIT   ,ICNT  ,IPP   ,   S1UM0011
     XJMAX  ,MW           ,PP    ,PTB   ,TRL1  ,TRL2  ,TRL3  ,TRL4  ,   S1UM0012
     XTRM1  ,TRM2  ,TRM3  ,TRM4  ,LRC   ,N3    ,FA    ,VGR              S1UM0013
C                                                                       S1UM0014
      TRM1(N)=0.                                                        S1UM0015
      TRM2(N)=0.                                                        S1UM0016
      TRM3(N)=0.                                                        S1UM0017
      TRM4(N)=0.                                                        S1UM0018
      DO 4 J=1,JMAX                                                     S1UM0019
      C1=A(J)*PTB(J+20,N)                                               S1UM0020
      TRM1(N)=TRM1(N)+C1                                                S1UM0021
      TRM2(N)=TRM2(N)+C1*PTB(J+30,N)                                    S1UM0022
      TRM3(N)=TRM3(N)+C1*PTB(J+40,N)                                    S1UM0023
 4    TRM4(N)=TRM4(N)+C1*((PTB(3,N)-PTB(J+30,N))**2+                    S1UM0024
     X     (PTB(4,N)-PTB(J+40,N))**2-PTB(J+30,N)/A(J)/PTB(7,N)*EFN(J,N))S1UM0025
      TRM4(N)=TRM4(N)*(GAMMA-1.)                                        S1UM0026
      RETURN                                                            S1UM0027
       END                                                              S1UM0028
$IBFTC S2UMP                                                            S2UM0000
      SUBROUTINE SUMP2(N)                                               S2UM0001
C                                                                       S2UM0002
      COMMON/NAMEP/P(50,100)                                            S2UM0003
      DIMENSION   A(10)    ,BS(10)   ,WPWG(10)           ,PP(50,100),   S2UM0004
     X PS(50,50) ,PTB(50,8),PW(150)  ,R(10)    ,TRL1(5)  ,TRL2(5)   ,   S2UM0005
     X TRL3(5)   ,TRL4(5)  ,TRM1(5)  ,TRM2(5)  ,TRM3(5)  ,TRM4(5)   ,   S2UM0006
     X FA(100)                                                          S2UM0007
      COMMON CAPN  ,CPL   ,CPS   ,CTOP  ,DL    ,DR    ,EK    ,EKW   ,   S2UM0008
     XEW    ,WPWG  ,GAMMA ,GMG0  ,HPL   ,HPS   ,IMAX  ,N1    ,N2    ,   S2UM0009
     XPS    ,PG0   ,PW    ,R     ,RCAP  ,RHOP0 ,SMP   ,TG0   ,TPM   ,   S2UM0010
     XA     ,BS    ,CC    ,CONST ,EOW1  ,EOW2  ,HIT   ,ICNT  ,IPP   ,   S2UM0011
     XJMAX  ,MW           ,PP    ,PTB   ,TRL1  ,TRL2  ,TRL3  ,TRL4  ,   S2UM0012
     XTRM1  ,TRM2  ,TRM3  ,TRM4  ,LRC   ,N3    ,FA    ,VGR              S2UM0013
C                                                                       S2UM0014
      TRM1(N)=0.                                                        S2UM0015
      TRM2(N)=0.                                                        S2UM0016
      TRM3(N)=0.                                                        S2UM0017
      TRM4(N)=0.                                                        S2UM0018
      DO 4 J=1,JMAX                                                     S2UM0019
      C1=A(J)*PTB(J+20,N)                                               S2UM0020
      TRL1(N)=TRM1(N)                                                   S2UM0021
      TRM1(N)=TRM1(N)+C1                                                S2UM0022
      TRL2(N)=TRM2(N)                                                   S2UM0023
      TRM2(N)=TRM2(N)+C1*PTB(J+30,N)                                    S2UM0024
      TRL3(N)=TRM3(N)                                                   S2UM0025
      TRM3(N)=TRM3(N)+C1*PTB(J+40,N)                                    S2UM0026
      TRL4(N)=TRM4(N)                                                   S2UM0027
 4    TRM4(N)=TRM4(N)+C1*((PTB(3,N)-PTB(J+30,N))**2+                    S2UM0028
     X     (PTB(4,N)-PTB(J+40,N))**2-PTB(J+30,N)/A(J)/PTB(7,N)*EFN(J,N))S2UM0029
      TRL4(N)=TRL4(N)*(GAMMA-1.)                                        S2UM0030
      TRM4(N)=TRM4(N)*(GAMMA-1.)                                        S2UM0031
      RETURN                                                            S2UM0032
       END                                                              S2UM0033
$IBFTC TAFN                                                             TAFN0000
      SUBROUTINE TAFN(X,Y,ANS)                                          TAFN0001
      ANS=TAN(.5*(ATAN(X)+ATAN(Y)))                                     TAFN0002
      RETURN                                                            TAFN0003
      END                                                               TAFN0004
$IBFTC WLPT                                                             WLPT0000
      SUBROUTINE WLPT                                                   WLPT0001
C     THE WALL POINT CALCULATION                                        WLPT0002
C                                                                       WLPT0003
        COMMON/ERR/IERR                                                 WLPT0004
      COMMON/NAMEP/P(50,100)                                            WLPT0005
      DIMENSION   A(10)    ,BS(10)   ,FJ(10)             ,PP(50,100),   WLPT0006
     X PS(50,50) ,PTB(50,8),PW(150)  ,R(10)    ,TRL1(5)  ,TRL2(5)   ,   WLPT0007
     X TRL3(5)   ,TRL4(5)  ,TRM1(5)  ,TRM2(5)  ,TRM3(5)  ,TRM4(5)       WLPT0008
     X,RT1(3)    ,RT2(3)   ,RT3(3)   ,VAL1(3)  ,CAPL(3)                 WLPT0009
     X,Q(3)                                                             WLPT0010
      COMMON CAPN  ,CPL   ,CPS   ,CTOP  ,DL    ,DR    ,INSERT,EKW   ,   WLPT0011
     XEW    ,FJ    ,GAMMA ,GMG0  ,HPL   ,HPS   ,IMAX  ,N1    ,N2    ,   WLPT0012
     XPS    ,PG0   ,PW    ,R     ,RCAP  ,RHOP0 ,SMP   ,TG0   ,TPM   ,   WLPT0013
     XA     ,BS    ,CC    ,CONST ,EOW1  ,EOW2  ,HIT   ,ICNT  ,IPP   ,   WLPT0014
     XJMAX  ,MW           ,PP    ,PTB   ,TRL1  ,TRL2  ,TRL3  ,TRL4  ,   WLPT0015
     XTRM1  ,TRM2  ,TRM3  ,TRM4                                         WLPT0016
C                                                                       WLPT0017
C     STEP 1 AVERAGEING FOR INITIAL CONDITIONS                          WLPT0018
C                                                                       WLPT0019
      DO 31 I=1,4                                                       WLPT0020
 31   PTB(I,3)=(PTB(I,5)+PTB(I,2))/2.                                   WLPT0021
C                                                                       WLPT0022
C     STEP 2 CALCULATION OF BARK                                        WLPT0023
C            BEGIN ITERATION LOOP                                       WLPT0024
C                                                                       WLPT0025
      DO 80 I=1,IMAX                                                    WLPT0026
      DO 21 N=2,3                                                       WLPT0027
      RT1(N)=SQRT (GAMMA*PTB(1,N))                                      WLPT0028
      RT2(N)=SQRT (PTB(2,N)*(PTB(3,N)**2 +PTB(4,N)**2 )-GAMMA*PTB(1,N)) WLPT0029
      RT3(N)=RT1(N)/RT2(N)                                              WLPT0030
 21   VAL1(N)=PTB(3,N)**2 *PTB(2,N)-GAMMA*PTB(1,N)                      WLPT0031
      DO 210 N=1,6                                                      WLPT0032
 210  PTB(N,8)=PTB(N,3)                                                 WLPT0033
 23   CAPK2=(PTB(2,2)*PTB(3,2)*PTB(4,2)+RT1(2)*RT2(2))/VAL1(2)          WLPT0034
      CAPK3=(PTB(2,3)*PTB(3,3)*PTB(4,3)+RT1(3)*RT2(3))/VAL1(3)          WLPT0035
      CALL TAFN(CAPK2,CAPK3,BARK)                                       WLPT0036
C                                                                       WLPT0037
C     STEP 3 CALCULATION OF    C, Z3, Q, R3, AND DECIDE PTS 6 AND 7     WLPT0038
C                                                                       WLPT0039
      C76=0.                                                            WLPT0040
 1    C76=C76+1.                                                        WLPT0041
      IF(C76-40.) 3,2,2                                                 WLPT0042
 2    WRITE (6,110)                                                     WLPT0043
 110  FORMAT(44H ERROR(10) UNABLE TO FIND PTS 6+7 ALONG WALL)           WLPT0044
      IERR=1                                                            WLPT0045
      RETURN                                                            WLPT0046
 3    PTB(5,6)=PW(MW-3)                                                 WLPT0047
      PTB(6,6)=PW(MW-2)                                                 WLPT0048
      PTB(5,7)=PW(MW-1)                                                 WLPT0049
      PTB(6,7)=PW(MW)                                                   WLPT0050
      C=(PTB(5,7)-PTB(5,6))/(PTB(6,7)-PTB(6,6))                         WLPT0051
      PTB(6,3)=(PTB(5,6)-PTB(5,2)+BARK*PTB(6,2)-C*PTB(6,6))/(BARK-C)    WLPT0052
      Q=(PTB(6,3)-PTB(6,6))/(PTB(6,7)-PTB(6,6))                         WLPT0053
      IF(Q-1.-EW) 7,7,4                                                 WLPT0054
 4    IF(PW(MW+2)) 5,6,5                                                WLPT0055
 5    MW=MW+2                                                           WLPT0056
      GO TO 1                                                           WLPT0057
 6    EOW1=1.                                                           WLPT0058
      GO TO 899                                                         WLPT0059
 7    IF(Q+EW) 8,9,9                                                    WLPT0060
 8    IF(MW-5) 81,81,82                                                 WLPT0061
 81   WRITE (6,111)                                                     WLPT0062
 111  FORMAT(42H0 ERROR(11) CHARACTERISTICS CROSS AT WALL.)             WLPT0063
      IERR=1                                                            WLPT0064
      RETURN                                                            WLPT0065
 82   MW=MW-2                                                           WLPT0066
      GO TO 1                                                           WLPT0067
 9    IF(PW(MW+2)) 12,10,12                                             WLPT0068
 10   IF(Q-1.+EW) 12,11,11                                              WLPT0069
 11   EOW2=1.                                                           WLPT0070
 12   EOW1=0.                                                           WLPT0071
      PTB(5,3)=PTB(5,6)+Q*(PTB(5,7)-PTB(5,6))                           WLPT0072
 899  IF(EOW1)881,400,881                                               WLPT0073
C                                                                       WLPT0074
C     STEP 4 CALCULATION OF GAS PROPERTIES P, RHO, AND U                WLPT0075
C                                                                       WLPT0076
 400  GMU2=PTB(2,2)*(PTB(4,2)*RT1(2)+PTB(3,2)*RT2(2))/(RT1(2)*VAL1(2))  WLPT0077
      GMU3=PTB(2,3)*(PTB(4,3)*RT1(3)+PTB(3,3)*RT2(3))/(RT1(3)*VAL1(3))  WLPT0078
      F2  =-RT3(2)*GAMMA*PTB(1,2)*PTB(4,2)/PTB(5,2)                     WLPT0079
      F3  =-RT3(3)*GAMMA*PTB(1,3)*PTB(4,3)/PTB(5,3)                     WLPT0080
      CAPJ2=PTB(1,2)+(PTB(6,3)-PTB(6,2))*(F3  *GMU3+F2  *GMU2)          WLPT0081
     X      /2.+(RT3(2)+RT3(3))*(                                       WLPT0082
     X      PTB(4,2)*(PTB(2,2)*PTB(3,2)+PTB(2,3)*PTB(3,3))-             WLPT0083
     X      PTB(3,2)*(PTB(2,2)*PTB(4,2)+PTB(2,3)*PTB(4,3)))/4.          WLPT0084
      CAPJ1=PTB(1,5)+(PTB(4,5)*(PTB(2,3)*PTB(4,3)+PTB(2,5)*             WLPT0085
     X    PTB(4,5))+PTB(3,5)*(PTB(2,3)*PTB(3,3)+PTB(2,5)*PTB(3,5)))/2.  WLPT0086
      G3=PTB(2,3)*PTB(3,3)                                              WLPT0087
      G5=PTB(2,5)*PTB(3,5)                                              WLPT0088
      H3=PTB(2,3)*PTB(4,3)                                              WLPT0089
      H5=PTB(2,5)*PTB(4,5)                                              WLPT0090
      CAPL(2)=RT3(2)*PTB(2,2)*PTB(4,2)                                  WLPT0091
      CAPL(3)=RT3(3)*PTB(2,3)*PTB(4,3)                                  WLPT0092
      Q(2)=  -RT3(2)*PTB(2,2)*PTB(3,2)                                  WLPT0093
      Q(3)=  -RT3(3)*PTB(2,3)*PTB(3,3)                                  WLPT0094
      S3=PTB(2,3)/(GAMMA*PTB(1,3))                                      WLPT0095
      S5=PTB(2,5)/(GAMMA*PTB(1,5))                                      WLPT0096
      PTB(3,3)=2.*(CAPJ1-CAPJ2)/(CAPL(3)+CAPL(2)+G5+G3+                 WLPT0097
     X         C*(Q(3)+Q(2)+H5+H3))                                     WLPT0098
      PTB(4,3)=PTB(3,3)*C                                               WLPT0099
      PTB(1,3)=CAPJ1-(PTB(3,3)*(G5+G3)+PTB(4,3)*(H5+H3))/2.             WLPT0100
      PTB(2,3)=PTB(2,5)+(S3+S5)*(PTB(1,3)-PTB(1,5))/2.                  WLPT0101
      NO=0                                                              WLPT0102
      DO 791 N=1,6                                                      WLPT0103
      CALL CRIT(PTB(N,3),PTB(N,8),NO)                                   WLPT0104
 791  NO=NO+NO                                                          WLPT0105
      IF(NO) 80,881,80                                                  WLPT0106
 80   CONTINUE                                                          WLPT0107
      I=IMAX                                                            WLPT0108
 881  ICNT=I                                                            WLPT0109
      RETURN                                                            WLPT0110
       END                                                              WLPT0111
